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PROJECT  SUMMARY 


A  self-c<Mitained  mesoscale  numerical  weather  prediction  system  which  can  generate  real¬ 
time  or  historical  simulations  on  a  high  performance  moderate-cost  workstation  coii^)uter  was 
developed  and  tested  in  this  project.  The  system  is  designed  to:  (1)  ingest  a  variety  of 
atmospheric  data  types;  (2)  combine  the  diverse  mixture  of  ingested  data  into  a  dynamically 
consistent  initialization  dataset  for  the  numerical  model;  (3)  generate  a  mesoscale  numerical 
simulation;  (4)  display  tlw  ouq)ut  tom  the  simulations  in  a  variety  of  graphical  and  tabular 
ftxmats;  and  (5)  be  ea^y  used  a^  reconfigured  through  a  graphical  user  interface. 

The  project  consisted  of  three  segments.  In  the  first  portion  of  the  project,  the 
wOTkstation-bas^  mesoscale  atmospheric  simulation  system  was  created  by  porting  a  version  of 
the  Mesoscale  Atmospheric  Simulation  System  (MASS)  tom  a  supercomputer  to  a  Stardent  750 
vector  i»ocessin^  workstation,  and  then  upgrading  the  system  by  implementing:  (1)  a  new  lateral 
boundi^  condition  scheme;  (2)  a  positive  definite  advection  formulation;  (3)  a  more 
sophisticated  surface  energy  and  moisture  budget  formulation;  (4)  a  four-dimensional  data 
assimilaticHi  system  based  on  a  Newtonian  relaxation  scheme;  and  (5)  a  method  to  enhance  the 
initialization  of  relative  humidity  through  the  use  of  satellite  image  data,  manually  digitized  radar 
(MDR)  reports,  pilot  reports  and  surface  cloud  observations.  In  the  second  portion  of  the  project, 
a  graphictti  user  interface  (GUI)  was  developed  to  permit  the  user  to  easily  reconfigure  the 
system  and  execute  simulations,  and  the  model  software  was  modified  to  increase  its 
computational  performance  on  the  workstation  computer.  The  final  segment  of  the  project 
consisted  of  the  execution  of  a  observation  system  simulation  experiment  (OSSE)  to  test  the 
performance  of  the  data  assimilation  system,  and  the  execution  of  36  real  data  simulations  to 
evaluate  the  performance  of  the  simulation  system  in  a  variety  of  environments. 

The  36-case  evaluation  experiment  indicated  that  the  modeling  system  was  able  to 
consistently  simulate  some  classes  of  phenomena  quite  well.  These  included  (1)  synoptic  or 


meso-a  scale  features  which  are  well-resolved  by  the  rawinsonde  network,  (2)  convective  events 
which  are  strongly  forced  by  larger  scale  circulations,  (3)  mesoscale  circulations  which  are  tied 
to  surface  features  that  are  well-represented,  e.g.  land/sea  breezes,  terrain-induced  circulations; 
(4)  surface  and  boundary  layer  flows  occurring  under  quiescent  (clear,  non-convective) 
conditions.  However,  the  evaluation  experiment  also  suggested  that  there  was  a  set  of 
phenomena  that  was  difficult  to  simulate  well.  These  phenomena  included  (1)  convective  events 
which  are  only  weakly  forced  by  larger  scale  circulations  (i.e.  quasi-barotropic  systems);  (2) 
features  which  result  tom  convective-scale  feedback  to  the  grid  scale  (e.g.  from  latent  heat 
release,  downdrafts,  etc.);  (3)  iiKSOscale  circulations  which  are  driven  by  more  subtle  surface 
gradients,  especially  soil  moisture;  (4)  circulations  driven  by  the  effects  of  cloud  boundary  on  the 
surface  energy  budget;  and  (S)  features  in  the  vicinity  of  large  terrain  gradients. 

In  addition  to  the  generation  of  real-time  mesoscale  simulations  for  military  applications, 
it  is  envisioned  that  the  workstation-based  mesoscale  simulation  system  could  have  a 
considerable  number  of  applications  in  the  private  sector  including  the  generation  of  tailored 
local  forecasts  for  business  and  government  operations  that  are  weatiier  sensitive.  A  number  of 
organizations  in  the  transportation,  recreation,  construction  and  utility  industries  could  benefit 
tom  high  resolution  short-term  weather  forecasts  that  could  be  provid^  by  such  a  system. 

Three  signiHcant  factors  should  permit  the  computational  and  meteorological 
performance  of  the  simulation  system  to  dramatically  increase  during  the  next  few  years:  (1)  the 
computational  power  of  workstation  computers  will  increase  and  thereby  permit  higher: 
resolution  models  with  more  complex  physics  to  executed  on  workstation-class  machines;  (2) 
model  initialization  data  with  mesoscale  and  cloud-scale  resolution  will  become  available  from- 
new  atmospheric  observing  systems;  and  (3)  the  physics  of  the  model  will  be  improved  as  the 
knowledge  of  mesoscale  and  cloud-scale  phenomena  increases.  Avaiiahiiitv 
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1.  Introduction 

Atmospheric  conditions  can  have  a  significant  impact  on  noilitary  operations  within  a 
battlefield  environinenL  Some  of  the  environmental  factors  which  can  impt^  tin  success  of  a 
milit^  operation  include:  (1)  the  state  of  the  ground  surface;  (2)  vertical  and  horizontal 
visibility;  (3)  the  temperature  of  the  ground  and  the  surface  air,  (4)  the  vertical  profile  of 
temperature  and  moisture;  and  (5)  the  vertical  profile  of  the  wind  direction  and  speed. 
Significant  precipitation  can  cause  the  soil  to  become  saturated  and  turn  to  mud  which  can 
adversely  imp^t  the  abili^  of  a  vehicle  to  quickly  traverse  an  area.  Clouds  and  fog  can  impose 
severe  restrictions  on  visibility  which  can  inhibit  aircraft  operations  and  even  ground-based 
activities.  Extreme  cold  or  heat  at  the  surface  of  the  earth  can  cause  problems  with  the  operation 
of  some  types  of  equipment  and  with  the  ability  of  personnel  to  efficiently  i^oim  critical  tasks. 
The  vertical  profiles  of  teinperature  and  moisture  can  impact  the  propagation  characteristics  of 
electromagnetic  and  acoustic  sensing  systems  and  can  prevent  or  enhance  the  ability  of  personnel 
to  acquire  critical  information.  The  dii^on  and  speed  of  the  wind  can  play  an  iiiq>ortant  role  in 
determining  the  marmer  in  which  substances  (some  which  may  be  toxic)  injected  into  the  air  as  a 
planned  or  inadvertent  result  of  explosive  detonations  are  dispersed.  The  nuuiner  in  which  such 
substances  are  dispersed  by  the  atmospheric  circulation  may  be  critical  to  the  survival  of 
personnel  in  an  area. 

Due  to  these  potential  effects  of  weather  on  a  variety  of  military  operations,  weather 
information  can  play  an  importwt  role  in  the  decision-making  process.  In  recent  years,  weather 
infonnation  for  military  o^rations  has  been  extracted  by  forecasters  from  the  output  of  large 
scale  numerical  weather  prediction  nKxlels  that  are  executed  at  centralized  military  forecast 
centers  such  as  the  Air  Force  Global  Weather  Central  (AFGWC)  in  Omaha,  Nebraska,  or  the 
Fleet  Numerical  Oceanographic  Center  (R^OC)  in  Monterey,  California.  These  forecasts  have 
been  quite  successful  at  providing  information  about  the  evolution  of  general  conditions  over  an 
area  of  interest.  However,  it  has  long  been  recognized  that  weather  conditions  can  vary 
significantly  on  short  space  and  time  scales  within  an  area  of  operations,  which  is  typically  a  few 
hundred  kilometers  on  each  side.  For  example,  heavy  rain  from  thunderstorms  can  create  very 
moist  soil  conditions  that  drastically  reduce  Ae  trafficability  over  an  area  which  is  less  than  100 
km  on  a  side,  while  soil  conditions  in  other  nearby  areas  remain  relatively  dry  and  favorable  for 
the  movement  of  militapr  vehicles.  Similar  variations  can  occur  in  the  other  meteorological 
paramett^  that  have  a  direct  impact  on  military  operations.  Thus,  it  would  be  desirable  to  have 
forecast  information  about  the  mesoscale  (10  -  500  km)  and  cloud-scale  (<  10  km)  variability  of 
the  important  parameters. 

A  logical  api^ach  to  the  prediction  of  the  mesoscale  and  cloud-scale  variability  of 
weather-related  conditions  in  the  battlefield  environment  is  to  extend  the  dynamical  simulation 
model  technology  that  has  proven  so  successful  in  generating  global  and  continental  scale 
forecasts  to  the  mesoscale  and  cloud-scale.  Mesoscale  and  cloud-scale  models  which  could  be 
used  for  this  purpose  have  been  under  development  by  several  universities  (c.g.  Perkey,  1976; 
Warner  and  Seaman,  1990;  Tripoli  and  Cotton,  1989)  and  federal  research  organizations  (Hodur, 
1987;  Mesinger  et  al.,  1988),  as  well  as  a  few  private  companies  Kaplan  et  1982c)  since  the 
late  1970's.  Research  investigations  have  shown  that  these  models  have  considerable  skill  in 
simulating  the  mesoscale  and  cloud-scale  variability  associated  with  some  types  of  atmospheric 
phenomena  (e.g.  Kocin  et  al.,  1985;  Zhang  and  Fritsch,  1988),  but  have  considerable  difficulty  in 
correctly  simulating  other  classes  of  mesoscale  or  cloud-scale  phenomena.  Nevertheless,  the 
perception  among  many  meteorologists  is  that  these  models  offer  the  potential  to  significantly 
improve  the  ability  to  produce  accurate  shon-term  forecasts  over  small  regions  in  at  least  some 
situations.  However,  these  models  must  be  executed  at  a  very  high  computational  rate  if  their 
output  is  to  be  available  in  real  time.  Until  recently,  it  was  only  possible  to  generate  forecasts  (as 
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opposed  to  historical  simulations)  from  these  models  by  using  a  high  performance 
supercomputer,  such  as  a  Cray  XMP  or  Cray  2,  which  cost  tens  of  millions  of  dollars.  In  fact, 
the  generation  of  real-time  cloud-scale  simulations  is  still  slightly  beyond  the  capability  of  the 
fastest  superconi^uters.  The  requirement  for  extremely  high  computational  power  meant  that  the 
cost  of  generating  real-time  mesoscale  numerical  simulations  was  very  high  and  that  the 
forecasts  could  be  generated  only  at  (^)erational  centers  which  possessed  such  computer  systems. 
In  most  instances,  this  was  not  feasible,  because  the  supercomputer  resources  at  these  centers 
were  already  totally  committed  to  the  ingestion  and  processing  of  the  real-time  atmospheric  data 
and  the  execution  of  global  and  regional  forecast  models.  Thus,  computational  resources  have 
^nerally  not  been  available  at  operational  centers  to  execute  real-time  high  resolution  mesoscale 
simulations  over  limited  areas. 

Recent  advances  in  computer  technology  have  provided  ah  altemadve  approach.  High 
performance,  moderate-cost  computer  workstations  now  have  the  computational  power  to 
execute  a  high  resolution  (e.g.  10  Im  grid  increment)  mesoscale  simulation  over  a  SOO  x  500  km 
area  in  time  to  have  a  significant  peri^  of  forecast  utility.  This  ciq^bility  raises  Ae  possibility 
of  executing  mesoscale  simulations  over  limited  areas  on  workstation  computers  at  local  sites 
where  forecast  information  is  needed.  Such  a  system  would  have  several  advantages  for  use  in 
military  applications.  First,  the  components  of  tiie  modeling  systems  could  be  selected  for  each 
application  based  upon  which  configuration  of  the  model  is  likely  to  perform  best  in  a  particular 
environment.  This  is  an  issue  because  certain  components  of  the  model  physics  are  not 
understood  sufficiently  well  to  formulate  a  general  mathematical  treatment  that  is  valid  for  all 
envircmments.  An  example  of  this  problem  is  the  parameterization  of  the  effects  of  sub-grid 
scale  moist  convection.  At  present,  there  is  no  universal  parameterization  that  works  equally 
well  for  all  grid  scales  and  all  environments.  Thus,  there  can  be  an  advantage  to  selecting  an 
appropriate  scheme  for  a  particular  application.  This  would  be  very  difficult  to  do  at  a  forecast 
cen^  which  must  execute  a  simulation  over  a  large  domain  which  encompasses  many  types  of 
environments.  However,  it  is  quite  feasible  for  loc^  area  simulations  on  a  computer  workstation. 
This  also  has  the  possibility  of  accelerating  the  computational  rate  of  a  simulation  if  some 
physics  can  be  omitted  from  a  simulation  in  a  particular  situation  without  significantly  degrading 
the  forecast  A  s^ond  advantage  of  the  local  workstation  system  is  that  the  system  could  be 
easily  adapted  to  ingest  local  datasets  that  might  not  be  available  to  the  operational  forecast 
center.  For  example,  the  military  has  (and  will  enhance)  the  capability  to  gatiier  meteorological 
data  throughout  die  battiefield  environment  through  the  use  of  surface  and  airlx>me  sensors 
which  accompany  operational  units.  This  data  could  be  transmitted  back  to  the  local  workstation 
computer  for  assimilation  into  the  local  simulations.  A  third  point  is  the  fact  that  the  local 
system  could  continue  to  ingest  local  data  and  generate  updated  simulations  even  if 
communications  with  the  central  forecast  center  were  lost  during  an  operation.  Finally,  the  local 
simulation  system  provides  a  large  degree  of  flexibility  to  the  decision-makers  which  require 
timely  and  easily-inteipreted  display  of  meteorological  information  for  use  as  a  Tactical  Decision 
Aid  (TDA).  Simulations  could  be  generated  by  the  workstation  system  on  demand  with  the 
latest  information.  Thus,  if  other  foctors  dictate  a  requirement  for  a  re-evaluation  of  the 
opmtional  plan,  updated  weather  forecast  information  could  be  generated  and  displayed  in  an 
easily  inter^eted  format  for  quick  access  during  the  decision-making  process.  Such  an  on- 
demand  service  would  be  difficult  for  a  central  forecast  center  to  satisfy,  since  it  generally  has  a 
commitment  to  provide  numerical  forecast  products  to  a  much  larger  area.  A  redirection  of 
resources  to  the  localized  simulation  could  compromise  the  quality  of  weather  information  for 
other  regions.  All  of  these  factors  suggest  that  it  would  be  desirable  to  construct  a  self-contained 
workstation-based  mesoscale  atmospheric  simulation  system.  This  was  the  basic  objective  of 
this  Small  Business  Innovation  Research  (SBIR)  project. 

It  was  envisioned  that  a  workstation-based  mesoscale  simulation  system  would  also  have 
a  considerable  number  of  applications  in  the  private  sector.  One  application  with  considerable 
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potential  is  the  generation  of  tailored  local  forecasts  for  business  and  government  operations  that 
are  weather  sensitive.  A  number  of  organizations  in  the  transportatimi,  recreatitm,  construction 
and  utility  industries  could  benefit  from  high  resolution  short-term  weather  forecasts  that  could 
be  provided  by  such  a  system.  Universities  could  also  utilize  the  woiitstation-based  simulation 
system  as  an  instructional  or  research  tool. 

The  fundamental  approach  to  the  creation  of  a  high  quality  workstation-based  simulation 
system  in  this  project  was  to  pmt  an  existing  mesoscale  atnnospheric  simulation  nradel  to  a  high 
performance  computer  workstation  and  to  enhance  the  system  by:  (1)  broadening  its  ability  to 
assimilate  a  diverse  mixture  of  atmospheric  data  types;  (2)  improving  the  physics  of  the  m^el 
by  incorporating  state-of-the-art  formulations  for  certain  key  components  of  the  atiiK>spheric 
physics;  (3)  selecting  an  appropriate  high  performance  moderate-cost  workstation  computer 
system  to  serve  as  the  computational  platform  for  the  simulation  system  and  then  optimizing  the 
inodel  to  execute  efficiently  on  this  system ;  and  (4)  developing  a  ^phical  user  interface  (GUI) 
to  permit  the  user  to  easily  and  quickly  configure  the  data  assimilation  and  simulation  model  for 
a  particular  application.  The  meeting  system  that  was  selected  to  port  to  the  workstation  was 
the  Mesoscale  Atmospheric  Simulation  System  (MASS)  model  (Kaplan  et  al.,  1982c).  The 
history  and  structure  of  this  modeling  system  is  described  in  Chapter  2  of  this  report. 

The  technical  feasibility  of  the  workstation-based  mesoscale  atmospheric  simulation 
system  was  demonstrated  during  the  Phase  I  portion  of  this  project,  by  executing  a  MASS 
simulation  over  the  Honduras/^icaragua  border  region  of  Ontral  America  on  an  Alliant  FX-1 
computer  system.  The  model  was  able  to  execute  at  a  sufficient  rate  to  permit  real-time  forecasts 
to  be  generated  over  small  domains.  Furthermore,  the  numerical  model  was  able  to  replicate 
realistic  meteorological  features  in  the  flow  field,  including  the  development  and  inland 
propagation  of  a  sea  breeze  convergence  front  and  intricate  mountain/valley  wind  circulations  in 
regions  of  complex  terrain.  Bursts  of  a  low-level  ’’dust"  tracer  and  plumes  of  an  upper-level 
"smoke"  tracer  were  injected  into  the  model  grid  domain  and  their  transport  was  successfully 
simulated.  The  advection  of  these  optical  obscurants  conformed  to  the  detailed  flow  patterns 
generated  by  the  locally-forced  atmospheric  circulations. 

The  objective  of  the  Phase  II  project  was  to  transform  the  modeling  system  into  a  tool 
that  could  efficiently  be  used  to  generate  quality  real-time  or  research  simulations  on  a  high 
performance  workstation  computer.  The  first  task  during  Phase  II  was  to  select  an  appropriate 
computational  platform  for  the  simulation  system.  The  Stardent  7S0  workstation  was  chosen 
because  it  had  a  custom  vector  processor  that  was  similar  to  those  found  on  supercomputers  such 
as  the  Cray  systems.  Since  the  MASS  code  was  originally  designed  to  execute  efficiently  in  a 
vector  processing  environment,  it  was  anticipated  that  it  would  perform  quite  well  on  the 
Stardent  system,  even  though  the  system's  scalar  processing  speed  was  not  significantly  faster 
than  other  workstation  computers  available  at  the  time  it  was  acquired.  Benchmark  tests  of  the 
MASS  code  proved  this  expectation  to  be  true.  Once  the  computational  platform  was  selected, 
four  types  of  enhancements  to  the  system  were  made;  (1)  the  numerical  algorithms  and  physics 
included  in  the  weather  prediction  model  were  improved  to  increase  the  accuracy  of  the 
mesoscale  simulations;  (2)  the  initialization  procedure  of  the  model  was  augmented  so  that  real¬ 
time  meteorological  data  from  local  battlefield  sensors  could  be  incorporated  into  the  forecasts; 
(3)  a  user-frienSy  graphical  user  interface  (GUI)  was  developed  in  order  to  allow  the  simulation 
system  to  be  operated  easily  with  minimal  training,  and  to  r^uce  the  probability  of  user  errors; 
and  (4)  the  simulation  system  software  was  optimized  so  that  the  time  to  execute  a  simulation  on 
a  workstation  computer  could  be  minimized.  This  report  describes  the  work  completed  in  each 
of  these  areas  during  Phase  II. 

An  overview  of  the  MASS  model  is  presented  in  Chapter  2.  The  improvements  to  the 
modeling  system  that  were  implemented  durinjg  this  project  are  described  in  Chapters  3, 4  and  5. 
The  improvements  included  the  implementation  of;  (1)  a  radiative  lateral  boundary  condition 
formulation  to  improve  the  simulation  of  the  flow  near  the  lateral  boundaries  of  the  model;  (2) 
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the  MPDATA  advection  scheme;  and  (3)  an  augmented  surface  energy  and  moisture  budget 
formulation.  The  implementation  of  a  four  dimensional  data  assimilation  capability  and  an 
enhanced  static  initiailization  capability  is  described  in  Chapter  6.  This  chapter  includes  a 
description  of  the  Newtonian  relaxation  (i.e.  nudging)  scheme  that  can  be  us^  to  assimilate 
rawinsonde,  surface  or  wind  profiler  data,  and  can  easily  be  extended  to  other  data  types.  An 
Observation  Simulation  System  Experiment  (OSSE)  was  conducted  to  analyze  the  impact  of 
different  data  types  on  the  data  assimilation  scheme.  The  design,  implementation  and  results 
from  this  experiment  are  described  in  Chapter  7.  The  modeling  system  was  evaluated  by 
executing  a  set  of  36  simulations  and  subjectively  and  objectively  verifying  the  results.  The 
results  of  the  evaluation  experiment  are  presented  in  Chapter  8.  The  development  and  design  of 
the  GUI  is  described  in  Chapter  9.  A  summary  of  the  wmk  completed  in  this  project  is  presentMl 
in  Chapter  10.  This  chapter  also  includes  a  discussion  of  the  areas  in  which  future  model 
develq)ment  is  required. 
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2.  Overview  of  the  MASS  Model 

The  ixMxleling  system  used  as  the  basis  for  the  workstatitm-based  mesoscale  simulation 
system  was  version  5  of  the  Mesoscale  Atmospheric  Simulation  System  (MASS).  The  Arst 
documented  version  (MASS  2.0)  of  this  model  was  described  by  Kaplan  et  al.  (1982c).  MASS 
was  originally  formulated  as  a  limited  area  hydrostatic  mesoscale  m^eling  system.  The  model 
finite-difference  equations  were  formulated  on  a  terrain-following  normalized  pressure 
coordinate  system  and  an  unstaggered  horizontal  grid  in  which  all  prognostic  variables  were 
computed  at  the  same  grid  point.  Sixth-order  accurate  horizontal  space  finite  differences  were 
used  to  calculate  spadal  derivatives  and  the  Euler-backward  time  marching  scheme  was 
employed.  A  30-case  sample  of  MASS  2.0  simulations  was  gatitered  during  the  spring  and 
summer  of  1982.  Tliis  sample  was  extensively  analyzed  by  scientists  at  the  Godard  Space 
Flight  Center.  The  synoptic  scale  forecast  skill  of  the  nnodel  was  evaluated  and  compared  to  the 
highest  resolution  m^el  then  operationally  available,  NMCs  LFM  model.  The  results  of  this 
study  were  reported  in  Koch  et  al.  (1985).  The  same  sample  of  cases  was  used  to  investigate  the 
ability  of  the  truxlel  to  predict  the  genesis  of  intense  mesoscale  convective  systems.  The  results 
of  that  analysis  was  reported  in  Koch  (1985). 

After  the  1982  test  period  the  MASS  2.0  model  was  upgraded  to  version  3.0.  The 
upgrades  were  described  in  Wong  et  al.  (1983).  The  improvements  included  additions  to  the 
surface  energy  and  moisture  budgets  and  the  implementation  the  Fritsch  and  Chappell  (1980) 
cumulus  parameterization  scheme.  A  version  of  the  Kuo-Anthes  (Anthes,  1977)  cumulus 
parameterization  schemes  was  implemented  in  1984.  The  6th-order  accurate  finite  difference 
approximations  were  replaced  with  4th-order  accurate  approximations  since  the  higher  order 
formulation  provided  a  minimal  benefit  at  a  significant  computational  cost.  The  MASS  3.0 
model  was  used  in  a  series  of  dynamical  case  studies  in  the  middle  1980’s.  The  MASS  model 
was  used  to  study  the  mesoscale  evolution  of  the  Grand  Island  tornado  case  of  June  3, 1980  by 
Kaplan  et  al.  (1982d),  Coats  et  al.  (1984)  and  Kaplan  et  al.  (1985).  Cbats  et  al.  (1984)  used 
MASS  to  investigate  the  effect  of  soil  moisture  gradients  on  the  evolution  of  the  mesoscale 
environment  in  this  event  21ack  et  al.  (1984)  stuped  the  effect  of  boundary  layer  fluxes  and 
deep  convective  processes  on  the  evolution  of  the  early  phases  of  the  east  coast  snowstorm  of 
February  10,  1983  (the  Megalopolitan  Snowstorm).  Uccellini  et  al.  (1983),  Uccellini  et  al. 
(1987)  and  Whitaker  et  al.  (1988)  used  the  Goddard  Space  Flight  Center  version  of  the  MASS 
model  to  investigate  the  role  of  jet  streak  dynamics  and  boundary  layer  fluxes  in  their  multi-year 
study  of  the  Presidents'  Day  snowstorm  of  February  18-19,  1979.  Kocin  et  al.  (1985) 
documented  the  performance  of  the  MASS  model  for  a  Washington  DC  snowburst  event  in 
March  of  1984.  Zack  and  Kaplan  (1987)  used  the  MASS  model  to  study  the  evolution  of  the 
severe  storm  environment  of  the  April  10, 1979  Wichita  Falls  tornado  outbreak  which  was  also 
the  first  field  day  of  the  AVE-SESAME  experiment. 

Version  4  of  the  MASS  model  was  formulated  in  1988.  The  Euler  backward  time 
marching  scheme  was  replaced  by  a  split  explicit  scheme  as  reported  by  Karyampudi  et  al. 
(1988).  In  addition  to  this  change,  a  version  of  the  Blackadar  boundary  layer  parameterization 
scheme  (Zhang  and  Anthes,  1982)  was  implemented.  2^ck  et  al.  (1988)  us^  this  version  of  the 
model  to  study  the  impact  of  synthetic  relative  humidity  data  derived  from  satellite  data  on  the 
short  term  simulation  of  convective  precipitation  over  the  Florida  peninsula.  Waight  et  al.  (1989) 
and  Waight  and  Zack  (1990)  used  MASS  to  simulate  the  evolution  of  convection  during  one  of 
the  cases  that  were  intensively  observed  during  the  Cooperative  Huntsville  Meteorological 
Experiment  (COHMEX).  Cram  et  al.  (1991)  utilized  version  4  of  MASS  to  conduct  an 
Observation  System  Simulation  Experiment  (OSSE)  to  test  a  scheme  to  retrieve  temperature  data 
from  wind  data  reported  by  the  experimental  network  of  four  wind  profiler  in  eastern  Colorado. 
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Version  S.O  of  the  MASS  model  was  developed  during  1991  and  1992.  This  version  of 
the  model  included  a  prognostic  grid  scale  moisture  scheme,  an  enhanced  surface  energy  budget, 
a  modified  Kuo  cumulus  parameterization  scheme  that  included  convective  scale  downdrafts  and 
a  nxne  comprehensive  long  and  short  wave  radiation  scheme. 

Ute  model  was  fu^er  developed  during  this  project  by  (1)  adding  a  four  dimensional 
dynamic  data  assimilation  capability  based  on  a  Newtonian  relaxation  scheme  and  the 
specification  of  vertical  profiles  of  latent  heating  and  moistening  from  cloud  and  precipitation 
observations.;  (2)  bro^ening  the  scope  of  the  surface  energy  and  moisture  budgets;  (3) 
implementing  a  radiative  lateral  boundary  condition  formulation;  (4)  incorporating  a  positive 
de^te  advection  scheme  for  selected  prognosdc  variables;  and  (S)  improving  the  computational 
efficiency  of  the  noodel  software.  At  the  end  of  the  project  the  version  of  MASS  which  included 
all  of  th^  enhancements  was  designated  as  version  S.5.  The  main  ctnnponents  of  MASS  5.5 
are  summarized  in  Table  2*1.  A  detailed  description  of  this  version  of  the  model  can  be  found  in 
the  MASS  Version  5.S  Reference  Manual  (MESO,  1993a).  The  components  that  are  shaded  in 
Tid>le  2*1  were  added  or  significantly  upgrad^  during  this  project  Copters  3  through  6  of  this 
report  provide  additional  detail  about  the  improvements  to  the  modeling  system  that  were 
completed  in  this  project 
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Table  2-1  Summary  of  MASS  5J. 


itialization  _ 


Autonuuic  calculation  of  tenain.  land/water  distribution,  land  use,  climatological  vegetation 
index  value  at  resolution  of  chosen  grid  domain. 

Accepts  several  types  of  NMC  gridded  data  as  flrst  guess  fields:  LPM,  NGM,  Gridded  Optimum 
Inteipoladon  (CJOI).  Data  from  previous  MASS  run  can  also  be  used  as  first  guess  fiekls. 

Soil  temperature  based  on  surface  temperature  averaged  over  the  previous  three  days. 

Re-analysis  using  either  Barnes  or  optimum  interpolation  objective  analysis  scheme  with 
significam  level  rawinsonde.  surface,  and  wind  profder  data. 

humidity  inofites 'decided  .-fiom  txafeot  doud 
-  mthmatly  dtgidzed  radar  (MDR)  data  and  mhared  and  vidbie 


umencal  Techniques 


3-D  prognostic  equations  for  u,  v,  T,  p,,  qy ,  q^  and 
Hydrostatic  assumption. 

Tenain  following  Op  (normalized  pressure)  vertical  coordinate. 

Arakawa  unsiaggered  "A'*  grid  on  a  stereographic  map  image  plane. 

Fourth-order  accurate  horizontal  space  differencing. 

Split-explicit  time  marching  scheme. 

Forward-backward  scheme  used  for  inertia-gravity  modes. 

Adams-Bashforth  scheme  used  for  the  adveclion  terms. 

'  Itositivnidefiiiheadyeeikmsdieroe  (MPDATA)av^  as  alternative  m  Adi^-Baidrfc^ 
Absorbing  upper  layer  can  be  used  to  damp  vertically  propagating  waves, 
f  l^masnb  dam  assuiuIatiQa  (mi^  wirul  pnrfiier,  sut&ce  and  radar  data. 


hysics 


Blackadar  high  resolution  PBL  parameterization. 

Detaikid  mtbtct  energy  and  mtdsture  but^ets  including  thiee-byer  sui&ce  liydRdogy  acheiire, 


Uses  high  resolution  (1  km)  USCS  land  use  and  vegetation  index  databases  to  determine  surface 
characteristics  such  as  roughness  height,  fraction  of  surface  covered  by  vegetation,  etc. 


•  Prognostic  equations  for  cloud  water  and  ice  (q^)  and  rain  water  and  snow  (q^)  (Diagnostic 

condensate  scheme  available  as  an  option). 

•  Simplified  parameterization  of  cloud  microphysicai  interactions. 

•  Kuo-type  cumulus  parameterization  with  moist  downdraft  physics  (Fritsch-Chappell  scheme 

available  as  an  option). 
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3.  Implementation  of  Radiative  Boundary  Conditions 

All  previous  versions  of  MASS  utilized  the  sponge  boundary  condition  originally 
proposed  by  Perkey  and  Kreiabecg  (1976).  In  this  formulation  the  tendencies  of  pro^ostic 
variables  from  an  external  source  are  blended  with  internally-generated  tendencies  within  a 
region  around  the  boundaries  of  the  domain. 

In  this  formulation,  the  tendency  of  a  prognostic  variable  at  a  point  near  the  boundary  is 
specified  by  the  relationship, 

^  =  W(I)^  +  [l-W(I)]^  (3-1) 

dt  dt  int  ^  ext 

where  the  subscript  "int"  denotes  an  internally-generated  model  tendency  and  the  subscript  "ext" 
denotes  a  tendency  from  an  external  source  such  as  a  larger  scale  model  or  an  analysis  of 
observational  data.  W  is  a  weighting  factor  which  determines  the  relative  contributions  of  the 
internal  and  external  tendencies  to  the  final  grid  point  tendency.  The  values  of  the  weighting 
coefficients  W(I)  are: 


W(I)  =1 


0.0  for  1  =  the  boundary  grid  points 
0.4  for  I  =  the  boundary- 1  grid  points 
0.7  for  1  =  the  boundary-2  grid  points 
0.9  for  I  =  the  boundary-3  grid  points 
1 .0  for  I  =  all  other  interior  grid  points 


I 


(3-2) 


If  these  weighting  coefficients  are  utilized,  then  the  boundary  point  is  completely  specified  by 
the  external  tendency  while  the  tendency  at  a  point  4A  from  the  boundary  is  completely 
determined  by  the  model  physics. 

Perkey  and  Kreitzberg  (1976)  demonstrated  that  the  effect  of  this  boundary  condition  is 
to  reduce  the  phase  velocity  of  a  disturbance  as  it  approaches  the  boundary.  This  transforms 
long-  and  medium-length  advective  and  gravity  waves  into  shon  waves  which  can  then  be 
removed  by  a  low  pass  filter,  thereby  giving  the  appearance  that  the  exiting  waves  simply  passed 
through  the  boundary.  This  formulation  performs  reasonably  well,  but  it  does  have  some 
drawbacks  for  mesoscale  simulations  on  workstation  computers. 

The  most  significant  negative  is  the  fact  that  at  least  4  points  adjacent  to  each  lateral 
boundary  must  be  used  to  apply  the  boundary  condition.  All  the  terms  in  the  equations  must  be 
solved  at  all  of  the  boundary  condition  points  with  the  exception  of  the  outermost  row  and 
column.  This  results  in  the  execution  of  a  sizable  number  of  computations  just  to  apply  the 
boundary  condition.  It  also  causes  the  size  of  the  interior  portion  of  the  domain  to  be  reduced 
since  a  total  of  8  points  in  each  horizontal  coordinate  direction  must  be  used  for  boundary 
conditions.  This  is  somewhat  acceptable  on  supercomputers  where  large  matrix  sizes  can  be 
employed.  However,  it  represents  a  considerable  fraction  of  the  domain  when  the  simulations 
must  executed  with  a  moderate  matrix  size  in  order  to  keep  the  processing  time  within  real-time 
constraints  on  a  workstation  computer.  For  example,  a  55  by  50  matrix  becomes  a  47  by  42 
matrix  of  interior  points. 

In  order  to  address  these  problems,  it  was  decided  to  implement  a  nonperiodic,  open 
lateral  boundary  condition  based  on  the  radiation  condition: 
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(3-3) 


at  '“dn 

where  X  is  a  prognostic  variable,  n  is  the  coordinate  perpendicular  to  the  boundary,  and  C  is  the 
phase  velocity  directed  normal  to  the  boundary.  C)rlanski  (1976)  developed  a  boundary 
condition  for  atmospheric  models  based  upon  this  condition.  In  the  Orlanski  scheme,  the  phase 
spee^  for  each  prognostic  variable  are  explicitly  calculated  at  the  nearest  interior  grid  point,  then 
applied  at  the  boundary  during  the  subsequent  timestep.  A  version  of  this  radiative  formulation, 
adapted  for  use  in  MESO's  TASS  (Terminal  Area  Simulation  System)  cloud  scale  nxxlel 
(doctor,  1985),  has  been  shown  to  allow  a  realistic  propagation  of  disturbances  through  the  open 
lateral  boundaries  with  a  minimum  of  wave  reflection. 


3.1  Orlanski  Radiative  Formulation 

In  the  Orlanski  (1976)  formulation,  the  values  of  the  dependent  variables  at  the 
boundaries  of  the  domain  are  determined  from  the  relationship 


^ 

•N  -V 


(3-4) 


where  x  is  any  prognostic  variable,  C  is  the  phase  velocity  of  waves  propagating  normal  to  the 
boundary  of  the  domain,  and  n  is  the  grid  coordinate  in  the  direction  normal  to  the  boundary.  In 
order  to  use  this  relationship,  the  phase  velocity  must  be  specified.  In  practice,  there  are  a 
number  of  wave  modes  which  are  in^inging  upon  the  boundary  at  a  particular  time,  each  with  its 
own  phase  velocity.  Consequently,  the  most  rigorous  application  of  (3-4)  requires  a 
decomposition  of  the  prognostic  variable  field  into  its  component  wave  modes  and  the 
calculation  of  a  phase  velocity  for  each  mode.  This  would  require  a  significant  amount  of 
computational  resources  if  it  had  to  be  done  at  every  model  timestep.  Fonunately,  a  simpler 
approach  has  proven  to  yield  reasonable  results  in  most  circumstances.  In  the  simpler  approach, 
a  composite  phase  velocity  is  estimated  by  solving  (3-4)  for  "C”  at  a  point  adjacent  to  the 
bound^  point: 


c  =  -  I-  .  ^ 

^  X‘  ~  X‘  At 

dn 

In  this  expression,  BP  denotes  the  boundary  point,  t  is  the  current  time.  At  is  the  length  of  a 

model  timestep,  and  An  is  the  grid  increment  in  the  direction  normal  to  the  boundary.  The  phase 
velocity  calculated  from  this  ex^ssion  represents  a  composite  value  which  does  not  correspond 
to  the  phase  velocity  of  a  particular  mode.  However,  it  should  be  close  to  the  mean  phase 
velocity  of  the  highest  amplitude  waves.  The  value  of  C  is  restricted  by  the  following 
requirements: 

(l)If-3z^i  /9z/9n  >  An/At  then 
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c  »  ^  . 

At 


(3^) 


This  repre^ts  the  maximuni  feasible  numerical  velocity.  That  is,  the  wave  would  move  one  grid 
increment  in  one  timestep. 

(2)IfO<-dx^  /dx/dn  <dn/Atthen 


C  =  (3-7) 

dn 

from  (3-4).  This  condition  represents  an  outward  propagation  at  less  than  the  maximum  speed. 
i3)lf  -dx/dt  Idxfia  <  0,  then  C  »  0  and 

X‘+^(BP)  =  At  ,  (3-8) 


where  {dpdt)nt  is  a  tendency  from  an  external  source  such  as  a  larger  scale  model  or  a  gridded 
observational  dataset.  This  case  represents  inward  propa'  .tion.  In  this  situation,  the  condition 
insures  that  no  interior  information  is  used  to  update  the  value  of  the  prognostic  variable  at  the 
boundary  point. 

32  Implementation  and  Testing 

Before  the  software  for  the  radiative  boundary  conditions  was  inserted  in  the  three- 
dimensional  version  of  the  MASS  model,  the  scheme  was  tested  in  a  simple  two-dimensional 
advection  program.  A  couple  of  two-dimensional  simulations  were  executed.  The  Hrst 
simulation  employed  a  simulation  grid  with  25  points  in  the  I  (left  to  right)  direction  and  20 
jmints  in  the  J  (up  and  down  direction).  The  ^d  spacing  was  specified  as  SO  km  and  the 
timestep  was  set  to  60  s.  A  tracer  value  of  10  units  was  ^)ecified  over  a  7  x  7  set  of  points 
centered  at  point  1=18,  J=10.  A  uniform  flow  of  25  m  s'^  along  the  I-^s  was  used  to  advect  the 
tracer  substance.  In  oirier  to  execute  the  simulations  quickly  on  a  Macintosh  fl  microcomputer, 
the  Ist-order  MPDATA  (donor  cell)  advection  scheme  was  used  for  the  advection  computations. 
This  resulted  in  a  more  diffusive  solution  than  would  be  the  case  had  the  3rd-order  F^DATA 
scheme  been  utilized  (see  section  5).  However,  this  should  not  have  had  a  significant  effect  on 
the  conclusions  about  the  performance  of  the  radiative  boundary  condition  software.  The  second 
simulation  was  identical  to  the  rirst  with  the  exception  that  the  I-dimension  of  the  domain  was 
expanded  from  25  points  to  35  points.  Thus,  the  boundary  column  (1=25)  on  the  first  run  was 
located  10  points  to  the  left  of  the  boundary  in  the  second  run.  Therefore,  the  tracer  values  for  the 
1=25  column  in  Run  1  (designated  the  “Boundary”  run)  were  determined  by  the  radiative 
boundary  conditions  while  in  Run  2  (termed  the  “Internal”  run),  they  were  the  result  of  the 
MPDATA  advection  calculations.  Both  simulations  were  execut^  for  a  total  of  720  tiroesteps 
which  is  equivalent  to  12  hours  of  simulated  tin^.  Figure  3- la  illustrates  the  history  of  the  tracer 
values  at  the  point  1=25,  J=10  for  both  Run  1  and  Run  2.  The  values  are  so  similar  that  only  one 
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curve  is  perceptible  in  the  figure,  even  though  data  from  both  runs  were  plotted.  The  diffoenc^ 
between  the  two  curves  is  shown  in  Figure  3- lb.  The  magnitude  of  the  difference  peaks 
s^pfoxiinately  30  to  40  timesums  before  the  tracer  value  peaks  at  the  boundary  point  At  25  m  S'‘ 
the  tracer  pe^  would  take  33.3  timesteps  to  move  across  one  SO  km  grid  cell.  Thus,  the  peak  in 
the  difference  values  occurs  at  the  time  when  the  peak  of  the  tracer  pool  is  ^proaching  the  point 
that  is  one  grid  interval  from  the  boundary.  An  explanation  for  this  behavitv  can  be  formulated 
by  recalling  that  the  phase  velocity  used  to  transport  the  tracer  from  the  interior  to  the  boundary 
point  is  determined  by 


C  =  ~  ,  (3-9) 

where  the  numerator  is  the  temporal  derivative  at  the  point  adjacent  to  oundary  and  the 
denominator  is  the  spatial  derivative  normal  to  the  boundary.  In  the  versiun  of  the  radiative 
boundary  condition  subroutine  that  was  implemented  in  the  MASS  code,  the  temporal  derivative 
is  determined  by  time  differences  of  the  tracer  values  one  point  in  from  the  boundary,  and  the 
spatial  gradient  is  calculated  from  the  difference  in  tracer  values  between  the  points  one  and  two 
gjid  intervals  in  from  the  boundary.  As  the  tracer  maximum  passes  through  the  area  bounded  by 
die  two  points  used  to  calculated  the  spatial  derivative  the  sign  of  the  derivative  must  change. 
Therefore,  the  spatial  derivative  history  curve  must  pass  through  zero.  The  software  is 
formulated  so  that  a  zero  value  can  never  be  used  in  the  calculation  of  the  phase  velocity. 
However,  v^  small  values  of  the  spatial  derivative  can  result  in  large  values  of  the  calculate 
phase  velocity.  The  only  limit  on  C  is  the  ratio  of  the  grid  spacing  to  the  timestep,  which  in  this 
case  is  833  m  s-^  Figure  3-lc  illustrates  the  variation  of  C  at  point  1=25,  J=10  for  the  Boundary 
run.  The  calculated  values  are  generally  very  close  to  the  actual  constant  advective  speed  of  25 
m  S'*  used  in  both  simulations.  However,  a  noticeable  perturbation  is  present  just  after  the  1 80th 
timestep,  when  the  tracer  peak  passes  between  points  1=23  and  1=24  and  the  calculated  value  of 
the  spatial  derivative  approaches  zero.  Despite  this  singularity  in  the  C  values.  Figure  3- la 
indicates  that  the  quality  of  the  simulation  at  ^e  boundary  point  is  still  quite  good. 

Another  view  of  the  evolution  of  the  Boundary  simulation  is  shown  in  Figure  3-2.  The 
initial  tracer  distribution  is  contoured  in  Figure  3-2a.  The  advecdon  of  the  tracer  downstream  is 
depicted  in  the  plots  for  120,  240  and  360  dmesteps  in  figures  3-2b,  3-2c  and  3-2d.  These  plots 
confirm  the  expectation  that  the  first  order  MPDATA  scheme  would  significantly  diffuse  the 
tracer  pool  during  the  course  of  the  simulation.  However,  it  can  also  be  seen  that  the  pool  passes 
through  the  bounda^  without  the  development  of  any  spurious  reflections  or  others 
perturbations.  In  addition  to  the  two  simulations  described  in  the  precnling  paragraphs,  several 
tests  were  run  with  different  magnitudes  for  the  advective  velocity  and  transpon  across  each  of 
the  four  sides  of  the  domain.  All  of  these  simulations  produced  satisfactory  results. 

The  Orlanski-type  radiative  lateral  boundary  condition  formulation  was  then 
implemented  into  the  three-dimensional  MASS  model.  An  experiment  was  conducted  to  verify 
that  the  boundary  condition  formulation  and  software  were  woricing  satisfactorily  and  to  compare 
the  perfOTinance  of  the  new  radiative  boundary  condition  to  the  Kreitzberg-Peikey  porous  sponge 
boundary  condition  that  was  previously  used  in  the  MASS  iiKxlel.  The  bound^  condition  test 
was  conducted  by  running  three  simulations.  Ail  three  simulations  were  initialized  with  an 
identical  analytical  dataset  generated  from  the  set  of  equations  formulated  by  Fritsch  et  al. 
(1980).  All  of  the  simulations  used  20  vertical  layers  extending  from  a  flat  surface  at  0  m 
elevation  to  100  mb  and  a  horizontal  grid  distance  of  20  km.  'Die  relative  humidity  was  set  to  a 
uniform  initial  value  of  1%  so  that  no  condensation  would  occur  during  the  simulation.  Each 
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simulation  was  executed  for  a  period  of  6  hours.  The  first  simulation  was  the  control  simulation. 
This  simulation  was  initialized  over  a  SO  by  40  horizontal  matrix  that  covered  the  entire  area 

shown  in  Hgure  3-3a.  The  second  simulation  was  integrated  over  a  40  x  30  matrix  that  covered 
the  area  depicted  by  the  inner  box  in  Figure  3-3a.  This  simulation  used  the  Kreitzberg-Perkey 
porous  sponge  lateral  boundary  condition  formulation  and  will  be  referred  to  as  the  “sponge” 
simulaticxi.  llie  lateral  boundaries  of  the  domain  for  the  second  simulation  are  witiiin  the  domain 
of  the  control  simulation.  Thus,  a  comparison  of  the  sponge  simulation  with  the  control 
simulation  provides  a  measure  of  the  impact  of  the  lateral  boundary  conditions  on  the  sponge 
simulation.  The  third  simulation  was  identical  to  the  second  simulation  with  the  exception  that 
the  radiative  lateral  boundary  condition  formulation  was  used.  This  simulation  will  be  referred  to 
as  the  “radiative”  simulation.  The  300  mb  height  and  wind  fields  at  the  end  of  each  6  hr 
simulation  are  shown  in  Figures  3-3b,  3*3c  and  3-3d.  These  fields  reveal  that  there  is  a  noticeable 
phase  error  in  the  sponge  simulation  which  is  much  less  significant  in  the  radiative  simulation. 
The  phase  error  is  revealed  by  the  position  of  the  main  ridge  axis  in  each  simulation.  In  the 
control  simulation  the  ridge  axis  extends  from  eastern  South  Dakota  southward  to  eastern 
Nebraska.  The  radiative  simulation  places  the  axis  in  almost  the  same  location.  However,  the 
sponge  simulation  places  the  ridge  axis  over  South  Dakota  well  to  the  west  of  its  position  in  the 
control  simulation.  An  examination  of  the  mass  and  momentum  fields  at  other  levels  indicates 
that,  in  general,  the  radiative  simulation  reproduced  the  features  of  the  control  simulation 
somewhat  better  than  the  sponge  simulation.  The  one  negative  aspect  of  the  radiative  simulation 
was  that  there  was  a  domain-scale  mass  depletion  during  the  course  of  the  simulation  that  was 
not  present  in  the  sponge  simulation.  One  manifestation  of  the  mass  depletion  is  that  the  average 
300  mb  heights  are  about  100  m  lower  in  the  radiative  simulation  than  they  are  in  the  sponge  or 
control  simulations.  The  mass  depletion  was  a  result  of  the  fact  that  the  radiative  toundary 
conctition  formulation  results  in  a  mixture  of  externally-specified  (inflow  points)  and  internally- 
specified  (outflow  points)  values  being  assigned  to  the  boundary  points  of  the  domain.  This  is 
not  a  problem  when  all  of  the  boundary  points  are  specified  from  a  single  external  source  as  is 
the  case  when  the  sponge  lateral  bound^  condition  is  employed.  This  problem  was  corrected  by 
enforcing  a  domain-scale  total  mass  constraint  as  part  of  the  radiative  boundary  condition 
formulation. 
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Figure  3-1  History  plots  from  the  two-dimensional  radiative  boundary  condition  test 
simulations:  (a)  tracer  values  from  the  point  (25,  10)  for  the  "Boundary"  and  "Internal” 
simulations.  Note  that  the  curves  for  the  two  smuUttions  are  virtually  coincident,  (b)  Difference 
in  tracer  values  at  point  (25, 10)  between  the  "Boundary"  and  "Internal"  simulations;  and  (c) 
calculate  phase  speed  "C"  from  the  point  (25, 10)  from  the  "Boundary"  run. 
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Figure  3-2  Contours  of  tracer  values  in  the  I-J  plane  from  the  "Boundary''  simulation  for:  (a) 
the  initial  time:  (b)  trfter  120  timesteps;  (c)  after  240  timesteps;  and  (d)  after  360  timesteps. 
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Figure  3-3  Results  from  a  lateral  boundary  condition  test  experiment  with  the  three-dimensional 
version  of  MASS:  (a)  depiction  of  the  small  (inner  box)  and  large  (outer  box)  domains  used  in 
the  experiment  and  300  mb  heights  (solid  lines,  m),  wind  vectors  and  isotachs  (dashed  lines, 
m  r^)  six  hours  after  the  initialization  time  for:  (b)  the  control  simulation;  (c)  the  simulation  with 
sponge  boundary  conditions;  and  (d)  the  simulation  with  radiative  boundary  conditions. 
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4.  Improvements  in  Surface  Physics 

During  the  course  of  this  project,  extensive  changes  to  the  naodel's  surface 
parameterization  schemes  were  made.  Si^ificant  changes  to  the  surface  energy  budget  scheme 
included  an  improved  evapotranspiration  scheme  and  a  better  treatment  of  the  subsoil 
temperature.  A  completely  revamped  longwave  and  shortwave  radiation  scheme  was 
implemented,  and  a  formulation  for  solar  radiaticm  on  sloping  terrain  surfaces  was  add^.  The 
planetary  boundary  layer  scheme  was  cleaned  up  and  recoded,  and  a  much  more  sc^histicated 
surface  hydroloi^  scheme  was  implemented.  The  improvements  in  several  of  these  areas  will  be 
discuss^  in  the  allowing  sections. 

4.1  Evapotranspiration  Scheme 

It  has  long  been  recognized  that  the  evaporation  over  land  surfaces  is  highly  variable  in 
both  space  and  time,  and  that  it  depends  heavily  on  local  surface  characteristics,  such  as 
vegetation  cover  and  structure,  which  are  not  well  understood  at  scales  larger  than  the 
microscale.  The  term  evapotranspiration  has  been  coined  to  describe  the  cmnbined  effects  of 
evapmation  from  various  surfaces  and  transpiration  from  plant  canopies.  Dr.  Joe  Russo  of  ZedX, 
Inc.  served  as  a  consultant  to  this  project  to  help  develop  a  comprehensive  evapotranspiration 
scheme  for  the  MASS  model.  There  are  separate  formulations  for  transpiration,  and  for 
evaporation  from  bare  soil  and  from  the  cover  reservoir,  which  represents  water  on  the  surface  of 
plants  and  other  surfaces  which  consists  of  intercepted  rainfall. 

The  latent  heat  flux,  E^rt,  represents  the  effects  of  evapotranspiration  from  three  separate 
sources: 

E|oi  “  ^soil  ”*■  ^cov  •  (^1) 

Eve,  the  transpiration  from  plants,  E^oy  is  the  evaporation  directly  from  the  top  layer  of  soil, 
and  Eeov  is  the  cover  layer  evaporation,  which  occurs  as  a  result  of  rainfall  interception.  An 
empiric^  formulation  for  Eyeg  has  been  developed  which  scales  the  actual  transpiration  to  the 
potential  evaporation  and  the  deep  layer  soil  moisture.  When  the  evaporative  demand  is  low, 
vegetation  can  transpire  freely  over  a  wide  range  of  soil  moisture  conditions.  As  the  evaporative 
demand  increases  however,  the  actual  transpiration  becomes  more  sensitive  to  soil  moisture 
conditions,  representing  in  a  simple  way  the  process  of  stomatal  closure  in  plants.  The 
transpiration  rate  is  parameterized  as: 

Eveg  =  min  j<l-®^vkvEp  ^^  2) 

Ev,g  is  not  allowed  to  exceed  the  incoming  solar  radiative  flux  density,  acknowledging  the  direct 
relationship  between  the  absorption  of  solar  radiation  and  transpiration  through  the  common 
mechanism  of  stomatal  control.  Ep  is  the  potential  evapotranspiration,  given  by 

Ep  =  p.  LvChV.  (qs  (Tg )  -  q.)  .  (4-3) 

where  Ly  is  the  latent  heat  of  vaporization,  q,  (Tg  )  is  the  saturation  vapor  mixing  ratio  at  the 
surface  temperature,  and  q,  is  the  atmospheric  vapor  mixing  ratio  near  the  surface.  The 
empirical  transpiration  coeffrcient  is  given  by 
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wheie  «i  and  M2  are 


ai  «  exp  (-.41457  +  2J527  E*  - 10.134  £•*  41.894  E*^) . 

M2  -  -27.914  +  80.100  E.  -  95.607  E.^  +  8.8831  E.^.  (4-5) 

E*  is  a  nocmalized  tranqmation. 


E.  «  .  (4-6) 

^R^iere  E _ is  a  transpiration  rate  which  is  taken  to  be  neariy  the  maximum  allowable  rate.  Ev» 

is  not  constrained  to  be  only  positive.  Figure  4-1  shows  the  variation  of  transpiration  with  b^ 
evaporative  demand  (potential  evapmxtion)  and  soil  moisture.  One  of  the  essential 
characteristics  of  the  scheme  is  that  it  mimics  the  tendency  of  plants  to  progressively  restrict 
evapontion  (by  stomatal  closure)  with  increasing  evaporative  demand,  sometimes  even  in  the 
presence  of  abundant  soil  moisture. 

There  are  several  serious  uncertainties  concerning  tranq)iiation.  First,  different  species 
of  plants  can  react  very  differently  to  similar  conditions.  Since  comprehensive  digital  data  on  the 

feographical  and  seasonal  distribution  of  species  is  hard  to  find,  this  is  a  serious  problem, 
econd.  many  of  the  translation  studies  which  have  been  conducted  have  been  performed  with 
plants  isolated  in  an  artificial  (experimental)  environment,  and  it  is  not  clear  how  to  exti^iolate 
the  transpiration  behavior  of  individual  plants  to  the  characteristics  of  a  complex  (peshaps  very 
heterogeneous)  plant  canopy. 

The  bare  soil  evaporation  is  a  modification  of  the  formulation  of  Noilhan  and  Pianton 

(1989): 


E|o3  —  (1  “■  ®v)  Pai-v  Qt  V*[hu<l»  (Tg  )  —  qj  ♦  (4-7) 

where  Gi,  is  the  fraction  of  the  non-vegetated  portion  of  the  model  grid  box  which  is  covered  by 
bare  soil,  which  is  considered  to  be  a  function  of  land  use  type.  The  relative  humidity  at  the 
ground  surface  is  given  by 


h„  = 


|0.5[l-cos(|^*  )],  ifw,<wfc  . 

I  1  ,  ifwiSwfc, 


(4-8) 


where  w^  is  the  field  capacity,  which  is  defined  as  0.75  w^. 

l^e  cover  evaporation  is  the  eviqioration  of  water  which  is  stored  in  a  cover  reservoir, 
made  up  of  snowfall,  dewfall  and  intercepted  rainfall.  As  will  be  discussed  in  the  hydrology 
section,  the  interception  parameterization  of  Mahfouf  and  Jacquemin  (1989)  is  used.  If  there  is 
intercepted  rainfall  present,  then  the  fraction  of  the  grid  covered  by  intercepted  rainfall  is 
assumed  to  be 


Ti  =  (; 


Wc 


WcmM  * 


(4-9) 
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where  Weis  die  iniooqMioa  reservoir  and  w —  is  the  auudmuin  imoant  d  wiser  allowed.  The 
cover  evapontioa  is  ih» 


Eeow 


min 


At 

CiC^Ep 


(4-10) 


where  the  top  expression  represents  the  evaporation  rate  which  would  occur  if  all  of  the 
inteicqited  mmstute  evi^xirated  in  one  surface  energy  budget  tiinesiqi.  and  the  lower  expression 
is  a  rate  which  is  diat  cover  moisture  can  ewqiorate  freely  at  the  potential 

rate  over  vegetated  areas.  In  the  presence  of  intetcqited  water,  transpiration  is  reduced  by  the 

factor  (1-Oi),  following  Mahfouf  and  Jacquemin  (1989).  If  the  water  viqior  gradient  is  reversed, 
then  condensation  (dew  formation)  is  allowed  to  occur  at  a  rate  of  OyEp,  with  the  constraint  that 
die  cover  reservoir  cannot  exceed  a  maximum  valiK. 


Figure  4-1  Transpiration  curves  as  a  function  of  soil  moisture  for  several  values  of  evt^orative 
demand.  All  transpiration  units  are  mm! hr.  The  soil  moisture  fraction  is  the  fractional  soil 
saturation,  not  the  volumetric  soil  moisture  fraction. 


4.2  Treatment  of  Subsoil  Temperature 

One  longstanding  model  initialization  problem  is  that  a  routine  dau  source  for  surface 
temperature  ("^n  temperature",  the  actual  temperature  of  the  ground  surface,  not  the  near- 
surface  atmospheric  temperature)  is  not  available.  The  traditional  solution  to  the  problem  has 
been  to  assume  that  the  surface  temperature  is  equal  to  the  atmospheric  temperature,  but  that 
assumption  is  clearly  incmrect  much  of  the  time.  To  test  the  sensitivity  of  the  scheme  to  initial 
surface  temperature,  a  series  of  runs  were  made  in  which  the  atnxispheric  temperature  remained 
constant,  and  the  surface  temperature  varied  from  255  to  280  K.  In  these  runs,  the  subsoil 
temperature  was  set  equal  to  5  K  less  than  the  surface  temperature.  Figure  4-2a  shows  the 
evolution  of  the  low  level  atmospheric  temperature  over  a  24  hr  period  for  two  initial  surface 
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temperatures.  There  is  a  dramatic  difference  in  the  nighttime  low  temperature  between  the  two 
runs,  even  larger  than  the  difference  between  the  initial  surface  temperatures.  This  suggests  that 
the  scheme  can  be  ve^  sensitive  to  the  initial  surface  temperature.  A  second  series  of  runs  with 
the  same  initial  variation  of  surface  ten:^)erature  but  with  the  subsoil  temperature  kept  constant  at 
262.5  K  is  shown  in  Figure  4-2b.  The  difference  in  the  nighttime  low  between  the  two  runs  has 
been  significantly  reduced,  and  there  is  no  significant  difference  in  daytime  highs.  Tlus  indicates 
that  it  is  the  variation  in  subsoil  temperatures  which  produces  most  of  the  difference  in  Figure  4- 
2a,  not  the  variation  of  surface  temperatures.  The  reason  for  this  is  that  the  ground  heat  flux  term 
in  the  surface  energy  budget  can  be  very  significant,  even  over  fairly  shon  periods.  A  method  of 
iiutializing  the  subsoil  temperature  from  an  average  atmospheric  temperature  over  the  previous 
several  days  was  developed  and  used  in  all  of  the  simulations  for  this  project 
The  ground  heat  flux  term,  H,„,  is  given  by 

^  -Tj)  ,  (4-11) 

where  x  is  the  number  of  seconds  in  one  day  and  T2  is  the  subsurface  “restoring”  temperature  in 
the  force-restore  method.  This  formulation  allows  T2  to  vary  with  the  surface  temperature,  on 
the  time  scale  of  about  one  day.  The  initialization  of  T2  is  not  obvious;  one  approach  is  to  define 
it  as  the  average  temperature  over  the  previous  few  days. 
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tlm*  (UTC) 


Fi^re  4^2  Diurnal  variation  of  low  level  temperature  for  two  one-dimensional  simulation 
initialized  with  the  same  analytical  vertical  temperature  profile,  surface  temperatures  of 260  and 
280  K,  and:  (a)  subsoil  temperatures  of  255  and  275  K,  respectively,  and  (b)  subsoil 
temperatures  of 262 5  K  for  both  simulations. 
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43  Solar  Radiation  on  Sloping  Surfaces 

The  effect  of  terrain  slope  was  added  to  the  calculation  of  incoming  solar  radiation.  The 
method  follows  the  discussion  in  Pielke  (1984).  The  solar  zenith  angle,  Z,  is  replaced  in  the 
incoming  shortwave  term  in  the  radiation  budget  with  a  modified  zenith  angle,  Z',  which  is  the 
angle  between  a  vector  normal  to  the  terrain  slope  and  the  solar  position. 

The  solar  zenith  angle  (Z)  is  de^ed  as 

cos  Z  =  cos  9  cos  5  cos  h  sin  5  sin  9  ,  (4-12) 

where  9  is  latitude,  5  is  the  solar  declination  angle,  and  h  is  the  hour  angle.  If  there  is  a  local 
terrain  slope,  then  the  following  calculations  are  made  to  incorporate  the  effect.  If  the  magnitude 
of  the  local  slope  is 


where  Zg  is  the  terrain  height.  The  azimuth  of  the  local  slope  is 


(4-13) 


(4-14) 


and  the  solar  azimuth  (which  varies  from  -Jt  to  +ji,  with  south  being  zero  and  west  of  south  being 
the  positive  direction)  is 


y,  =  sin-1  (gas§  anJi) . 
'  sin  Z  ' 

The  zenith  angle  corrected  for  slope  is  then 


(4-15) 


cos  Z'  *  cos  Oz  cos  Z  +  sin  Cz  sin  Z  cos  (Ys  -  Yz)  •  (4-1 6) 

Care  must  be  taken  that  no  division  by  zero  occurs  in  (4-14)  and  (4-15). 

To  show  the  effect  of  slope  on  incoming  solar  relation,  the  differences  between  the 

radiation  received  on  level  terrain  and  surface  sloping  (a  =  .05)  to  the  east  and  south  at  32*  N 
latitude  at  both  solstices  is  shown  in  Figure  4-3.  The  differences  are  large  enough  to  be 
significant  in  areas  of  mountainous  terrain,  especially  when  the  model  is  used  at  high  resolution 
(e.g.  10  km).  Although  this  formulation  is  a  pan  of  the  MASS  model,  no  three-dimensional 
simulations  with  the  terrain  slope  effects  were  performed,  because  a  high  resolution  terrain 
dataset  has  not  yet  been  incorporated  into  the  preprocessor.  At  the  end  of  the  project,  a  1  km 
terrain  dataset  which  would  allow  for  reasonable  slope  calculations  for  high  resolution  MASS 
grids  had  been  obtained,  and  a  relatively  modest  effon  is  needed  to  integrate  it  into  the 
preprocessor. 
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4.4  Hydrology  Scheme 

As  a  pan  of  the  effort  to  improve  the  model's  evapotranspiradon  scheme,  it  was  necessary 
to  construct  a  complete  hydrology  fnunework.  The  soil  hydrology  scheme  is  based  on  the  noodel 
of  Mahn  and  Pan  (1984).  Hie  soil  is  divided  into  two  layers  -  a  shallow  S  cm  layer  at  the 
surface  and  a  la^  fnmi  5  cm  to  30  cm  deep,  which  is  assumed  to  ctmtain  the  majority  of  plant 
roots.  An  additional  reservoir  of  moisture  is  parameterized,  a  "cover"  moisture  reservoir  which 
retains  intercepted  rainfall.  The  model  structure  and  processes  considered  are  depicted  in  Figure 
4-4.  Each  mo^l  surface  grid  box  has  been  assigned  a  land  use  type,  a  soil  type,  and  a  value  of 
fractional  vegetation. 

4.4.1  Soil  moisture  budgets 

With  the  notation  that  layer  1  is  the  shallow  layer,  and  layer  2  is  the  deep  layer,  soil 
moisture  budgets  for  the  two  layers  are 


—  =  :J-(INFIL-DIFF~CONi-E«,a)  . 

dt 

—  =  l--,(DIFF-CONi--CON2-Eveg)  , 
dt  (Z2-2l) 


(4-17) 


where  wi  and  W2  are  volumetric  soil  moisture  fractions  for  the  two  layers,  and  zi  and  Z2  are  the 
depths  of  the  bottom  of  each  layer  (S  cm  and  30  cm).  All  of  the  terms  in  parentheses  in  (4-17) 
have  units  of  m  s*>.  Many  of  the  expressions  for  the  terms  arc  from  McCumbcr  and  ihielke 
(1981). 

The  first  term  in  the  wi  equation  is  the  infiltration  of  precipitation  into  the  top  layer  of 
soil.  The  infiltration  is  the  precipitation  rate  minus  the  rainfall  interception: 

INFIL  =  PRECIP-ICEPT  .  (4-18) 

The  infiltration  rate  is  constrained  to  not  exceed  a  maximum  rate,  which  is  related  to  the 
properties  of  the  soil  and  the  degree  of  saturation  of  the  soil: 

INFIW  =  +  .  (4-19) 

T 

where 


Dj^ 


h  K^at 

'Vjai 


(4-20) 


and  b,  Kut  and  ysu  are  functions  of  soil  type. 

Diffusion  between  layers  assumes  that  the  soil  moisture  values  occur  at  the  midpoint  of 
the  layers,  with  a  linear  vertical  gradient: 
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4.4.2  Rainfall  interception  and  cover  moisture  reservoir 

The  parameterization  for  rainfall  interception  follows  that  of  Mahfouf  and  Jacquemin 
(1989).  The  prognostic  equation  for  the  cover  moisture  is 


=  ICEPT-Ecov  ,  (4-25) 

dt 

where  Wc  is  the  cover  moisture  with  units  of  length,  which  can  be  thought  of  as  the  depth  of 
moisture  on  the  surface.  The  interception  term  is  given  as 

ICEPT  =  max 

where  Oy  is  the  fractional  vegetation,  PRECIP  is  the  precipitation  rate,  and  where  Cinun  is  a 
minimum  interception  constant,  which  is  a  function  of  land  use  type.  Cjmin  is  nonzero  for  those 
land  use  types  which  are  assum^  to  intercept  rainfall,  even  when  Ae  NDVI  is  very  low,  such  as 
forests  and  rangeland.  For  example,  deciduous  forests  may  have  low  NDVI  in  the  winter,  but 
rainfall  interception  will  still  take  place  even  with  all  ^e  leaves  gone,  although  at  a  lower  rate. 
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The  ejmiessioii  for  the  cover  Uw  evapontioo  is  given  in  (4-10).  The  moisture  in  die  cover 
reservcSr  is  not  allowed  to  exceed  a  maximum  value,  parameteriaed  as 


OvCihoux,  ifOv>0 
0.01  Cl  hgixx ,  if  Ov  **  0  * 


(4-27) 


where  Ct  is  a  land  use  dependent  constant,  which  can  vary  from  0  to  about  3  for  forested  land. 
This  parameter  replaces  the  leaf  area  index  in  Mahfouf  and  Jacquemin.  hmu  is  a  maximum 
reservoir  value;  die  value  of  0.2  kg  m*^  from  Mahfouf  and  Jacquemin  is  used. 


4.4.3  Snow  cover 

A  parameterization  for  the  effects  of  snow  cover  frdlows  that  ci  Segal  et  al.  (1991).  The 
most  important  effect  is  a  strnig  increase  in  the  shonwave  albedo.  An  eiqiression  for  the  albedo 
is 


a  *  0.5  (a«o  +  otjo ) + 0.32  f(Z)  ,  (4-28) 

where  is  the  albedo  for  solar  radiation  with  a  wavelength  ^  0.7  mm  (0.95),  and  aio  is  the 
albedo  for  solar  radiation  with  wavelengths  >  0.7  mm.  The  last  term  represents  the  additional 
reflectance  of  snow  surfaces  at  large  soliff  zenith  angles; 


f(Z) 


0. 


i( _ t±J _ 

b'l  't-2bcosZ 
i/ _ t±J _ 

*»(l  +  2bcos80® 


.  Z^dO® 

,  60®<Z<80® 

.80®SZ 


(4-29) 


where  Z  is  the  zenith  angle,  and  b  is  a  ctmstant  set  to  2,  as  in  Se^,  et  al.  Several  other  surface 
variables  are  changed  for  snow  cover,  superseding  values  otherwise  derived  from  land  use  type. 

The  heat  cqiacity  is  set  to  9.305  x  10-^  K  m^  J-*,  a  value  for  fresh  snow  from  Oke  (1978).  The 

longwave  emissivity  is  set  to  0.99,  also  from  Oke.  The  roughness  length  is  set  to  1  x  1(H  over 
some  land  use  type  (agricultural  land,  rangeland,  and  barren  land),  but  is  unchanged  for  other 
types  (urban  lan^  forests,  open  water,  wetlands)  on  the  assumption  Aat  large  obstacles  dominate 
die  roughness  length.  With  snow  cover,  the  fractional  vegetation  is  set  to  zero,  so  that  any 
vegetative  effects  are  eliminated. 

To  test  the  snow  cover  paramterization,  a  similar  series  of  runs  were  made.  Figure  4-5 
demonstrates  that  snow  cover  effectively  decouples  the  atmosphere  from  the  soil  temperature  due 
to  the  insulatmg  effect  of  a  complete  snow  cover,  and  that  the  initial  surface  temperature  is 
unlmpt»tant,  since  an  equililvium  is  reached  very  quickly  due  to  the  low  heat  capacity  of  snow. 
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Figure  4-4  Schematic  diagram  of  the  MASS  model  hydrology  scheme. 
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tlm*  (UTC) 


Figww  4-5  DUtmal  variadon  cf  low  level  temperature  for  two  one-dimensional  simulation 
initudited  with  the  same  anaiydcal  vertical  wnperature  profile,  surface  temperatures  of 260  and 
280  K,  and  vridt  a  complete  snow  cover. 
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5.  Positive  Definite  Transport  Scheme 

When  simulating  nwist  processes,  it  is  imperative  to  resolve  strong  humidity  and  cloud 
water,  cloud  ice,  rain  water  and  snow  gradients.  Version  5  of  the  MASS  model  utilizes  a  ^lit- 
explicit  time  integration  scheme  (a  forward-backward  scheme  for  inertial-gravity  modes  and  the 
Adams-Bashforth  scheme  for  the  advective  modes)  and  fourth-order  accurate  advection  and 
diffusion  qierators.  Despite  these  measures,  the  forecast  fields  still  suffer  from  false  numoical 
dispersion  (Gibbs  oscillations)  and  produce  small  negative  amounts  of  water  vapor  and  cloud 
water,  cloud  ice,  rain  water  and  snow.  In  practice,  these  artificial,  high  frequency  ripples  are 
filtered  out  of  the  solution  by  the  lateral  ditnision  operator.  Any  residual  negative  values  of  the 
mixing  ratio  of  these  quantities  were  then  reset  to  a  small  ^sitive  amount  However,  this 
indiscriminately  smears  out  the  solution  so  that  the  fine  scale  structure  in  the  field  is  often  lost. 
Hiis  problem  was  adebessed  during  Phase  II  by  implementing  a  positive-definite  non-diffusive 
advection  scheme  into  the  simulation  model. 

5.1  Review  of  Candidate  Schemes 

Three  numerical  techniques,  each  developed  for  modeling  fluid  flows  characterized  by 
strong  shocks  and  discontinuities,  were  consider^  as  potential  candidates  for  implementation: 
(1)  Flux-Corrected  Transport  (FCT)  (Boris  and  Book,  1973;  Zalesak,  1979),  (2)  the 
Smolaikiewicz  Multidimensional  Positive  Definite  Advection  Transp(»t  Algorithm  (MPDATA) 
scheme  (Smolaridewicz,  1983a),  and  (3)  the  Piecewise  Parabolic  Method  (PPM)  (Carpenter  et 
al.,  1990).  Each  of  these  techniques  is  characterized  by  a  lack  of  numerical  dispersion  and  a  low 
level  of  inherent  diffusion.  To  illustrate  the  concepts  behind  these  three  numerical  methods, 
consider  the  continuity  equation  written  in  flux  form: 

$+V-(V<p)  =  0  (5-1) 


or 

^+V.F.O  (5-2) 

where  <p  is  the  material  substance  which  is  to  be  advected  through  the  model  domain  and  F  is  the 
flux  of  this  substance. 

The  flux-corrected  transpon  technique  separates  the  flux  divergence  calculation  into  two 
steps.  First,  a  time-advanced  solution  is  calculated  using  a  low-order  Effusive  scheme,  typically 
first-order  upstream  finite  differencing.  Then,  in  an  attempt  to  restore  the  shape  of  the  mass  field 
to  its  pre-diffused  form,  antidiffusive  fluxes  are  applied,  which  are  defined  to  be  the  difference 
between  the  high  (usually  fourth  or  sixth)-order  "ripple  producing"  fluxes  and  the  original 
diffusive  fluxes.  Adding  these  antidiffusive  fluxes  at  full  strength  is  equivalent  to  replacing  the 
low-order  fluxes  by  those  from  the  high-order  scheme.  If  such  a  prescription  creates  new  extrema 
(peaks  or  valleys)  in  the  solution,  then  the  antidiffusive  fluxes  are  systematically  reduced.  Thus, 
it  is  possible  to  re-concentrate  the  mass  at  each  grid  point  without  causing  overshooting. 
Because  FCT  virtually  eliminates  numerical  diffusion,  it  is  extremely  effective  in  advecting  even 
cusp-like  features  such  as  shocks  of  square  waves  without  a  serious  degradation  in  their  form. 
Mattocks  and  Bleck  (1986)  have  successfully  used  FCT  to  prevent  the  generation  of  negative 
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layer  thickiMSses  in  an  isentropic  channel  model  and  to  accurately  locate  the  intersection  of  the 
isentn^s  with  the  pound. 

Smolarldewicz  (1983a,  1983b)  developed  a  positive  definite  advection  scheme  which  has 
small  implicit  diffusion  and  a  lower  ccmiputational  cost  than  FCT.  The  scheme  was  further 
reflned  by  Smolarkiewicz  and  Clark  (1986)  and  was  given  the  name  of  Mulddimensional 
Positive  Definite  Advection  Tianspm  Algorithm  (MPDATA).  Smolarldewicz  retdized  that  the 
conventional  one-sided,  upstream-differenced  analog  of  the  continuity  equation  (written  in  flux 
form  for  one  spatial  dimension); 


9F 

dt’'“  dx 


where 


_  I  Uj+iiPi  ,  if  u>0 
2  \«i+l<Pi+l  .  tf  U<0 


(5-3) 


is  actually  a  centered  difference  representation  of  the  advection/diffusion  equation: 


(S-4) 


where  Kimpi  is  an  "implied"  diffusivity.  Smolarkiewicz  counteracts  this  implied  diffusivity, 
without  sacrificing  positive  definiteness,  by  following  each  advection  step  with  a  "negative 
diffusion"  or  conective  step: 


dtp  3 


3q>, 


(5-5) 


In  order  to  ensure  that  the  solution  remains  positive-definite,  this  can  be  re-formulated  as  an 
advection  equation: 


dtp  d  .  ^ 


(5-6) 


with  the  "antidiffusive  advection  velocity"  set  to  zero  when  the  transported  material  is 
completely  depleted  at  a  grid  point: 


Ud  = 


Kimpi  dy 
tp  dx 
0 


if  tp>0 
if  tp  =  0  ■ 


(5-7) 


Snnolarkiewicz  noted  that  the  restorative,  antidiffusive  effect  of  the  second  step  can  be  enhanced 
by  multiplying  the  advection  velocity  by  a  "correction  coefficient"  slightly  larger  than  1  or  by 
repeating  the  correction  step  during  each  iteration.  He  presented  impressive  results  from  solid- 
body  rotation  simulations  for  the  multidimensional  and  time-splitting  forms  of  the  equations.  The 
original  mass  field  retained  its  shape  while  sharp  gradients  were  maintained. 

The  third  numerical  technique  which  was  considered  for  implementation  into  the  MASS 
model  was  originally  developed  for  simulating  astrophysical  phenomena.  Known  as  the 
piecewise  parabolic  method  (PPM),  it  consists  of  a  rather  unconventional  approach;  dependent 
variables  are  represented  by  monotonic  parabolas  fit  to  each  grid  interval  instead  of  discrete  grid 
point  values.  Unlike  standard  curve-fitting  and  global  spectral  techniques,  each  parabola  is 
uniquely  constructed  for  a  specific  grid  box  and  each  grid  interface  is  considered  to  be  a 
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discontinuity.  Nonlinear  fluxes  across  grid  boundaries  can  therefcxe  be  calculated  explicitly  from 
the  analytic  (parabolic)  functions  using  Riemann  integration.  An  eloquent  description  of  the 
method  is  provided  by  Carpenter  et  al.  (1990).  The  authors  first  present  an  example  which  shows 
how  parabolas  are  constructed  for  an  analytic  transcendental  function  (see  Figure  S-1). 

Essentially,  this  procedure  requires  that  values  for  the  three  coefficients  (9o>  <pi>  92) 
determined  so  that  a  unique  quadratic  representation  of  the  dependent  model  variable: 

<p(x)  =  <po  +  <Pix+<P2x2  (5-8) 

can  be  defined  for  each  grid  box.  The  first  step  is  to  generate  lower-order  piecewise  linear 
functions  by  calculating  slopes  at  the  left  edge,  center,  and  right  edge  of  each  grid  interval,  then 
taking  the  slope  with  the  smallest  magnitude.  If  the  grid  box  or  "zone”  average  of  the  variable: 

Jrx-f  Ax/2 

(pdx  (5-9) 

x-Ax/2 

is  an  extremum,  then  the  slope  is  set  to  zero.  After  the  set  of  connected  lines  is  assembled,  the 
second  step  is  to  construct  a  first-guess  parabola  with  a  unique  cubic  curve  fit  from  the 
neighboring  zone  averages  and  slopes  which  surround  the  left  and  right  edges  of  the  grid  interval. 
The  first,  second,  and  third  derivatives  of  the  provisional  solution  are  then  calculated  to 
determine  whether  any  "contact  discontinuities”  (hydraulic  jumps)  exist  within  each  zone.  If 
some  are  found,  then  the  parabola  is  "steepened"  to  avoid  smearing  out  sharp  gradients.  Finally, 
any  undershooting  or  overshooting  is  eliminated  by  flattening  out  the  parabola,  if  necessary,  to 
prevent  the  creation  of  new  extrema  and  preserve  the  monotonicity  (positive  definiteness)  of  the 
solution. 

A  schematic  which  shows  how  the  PPM  advection  process  works  is  displayed  in  Figure 
5-2.  First,  zone  (grid  box)  averages  are  determined  by  integrating  the  initial  step  function 
distribution  over  Ae  width  of  each  grid  interval.  Parabolas  are  constructed  within  each  zone,  as 
previously  descriM,  then  the  entire  form  is  translated  to  the  right  during  the  advection.  Next,  the 
two  parabolas  which  now  lie  within  a  given  zone  are  integrated  analytic^ly  to  deteimine  the  new 
zone  averages  (step  function  distribution).  These  new  zone  averages  are  used  as  initial  data  for 
^e  subsequent  time  step,  in  which  a  new  piecewise  parabolic  representation  is  calculated.  Like 
its  progenitors,  the  piecevase  parabolic  method  is  positive  definite,  retains  the  integrity  of  steep 
gradients,  and  is  characterized  by  a  nearly  complete  lack  of  computational  diffusion.  In  addition, 
PPM  yields  solutions  of  similar  accuracy  at  half  the  spatial  resolution  when  compared  with 
conventional  finite  difference  methods. 
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THE  PEKIEWISE  PARABOLIC  METHOD 
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Figure  5-1.  A  graphical  illustration  which  shows  how  unique  parabolas  are  constructed  for  each 
grid  box  to  numerically  represent  an  analytical  solution  in  the  piecewise  parabolic  method 
(PPM).  In  this  example  the  analytical  solution  is  a  hyperbolic  tangeru  function.  It  is  denoted  by 
the  heavy  black  line  in  the  first  panel  of  the  figure  (erfter  Carpenter  et  al.,  1990). 
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PIECEWISE  PARABOLIC  METHOD  (PPM) 
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Figure  5-2.  A  schematic  which  illustrates  how  the  piecewise  parabolic  method( P PM) 
process  works.  The  initial  analytical  function  is  denoted  by  the  heavy  black  line  in  the  first  panel 
of  the  figure  (cfter  Carpenter  et  al..  1990). 
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During  the  Phase  II  project,  a  survey  of  the  scientific  literature  revealed  that  PPM  was  too 
computationally  expensive  and  too  time-consuming  to  code  from  scratch.  Experience  with  the 
flux-corrected  transport  method  (Mattocks  and  Bleck,  1986),  a  second  possible  choice,  had 
shown  it  to  be  extremely  effective  in  adveedng  shock  waves  or  other  contact  discontinuities  with 
minimal  degradation  in  their  form.  However,  FCT  sometimes  generates  artificially  steep 
gradients  in  regions  where  none  should  exist,  typically  at  the  leadmg  edge  of  the  waveform. 
Other  numerical  modelers,  most  notably  Smolaridewicz  (1983a,  1983b),  are  critical  of  FCT 
because  it  can  become  computationally  expensive,  depending  on  how  well  the  flux-limiting 
decision  process  is  formulated.  Therefme,  MESO  selected  the  ^^DATA  scheme  as  the  positive - 
definite  advection  scheme  to  be  implemented  in  the  MASS  model. 

5.2  One  and  Two  Dimensional  Tests  ofMPDATA 

Befme  the  MPDATA  scheme  was  implemented  into  the  three-dimensional  version  of  the 
MASS  model,  a  number  of  one-dimensional  and  two-dimensional  experiments  were  executed  to 
test  and  evaluate  the  scheme  under  a  simple  set  of  conditions.  The  objective  of  these 
experiments  was  to:  (1)  acquire  a  comprehensive  understanding  of  the  scheme's  performance,  (2) 
determine  the  computational  performance  of  each  of  the  scheme's  several  options,  and  (3)  ensure 
that  the  computer  code  us^  to  implement  the  scheme  was  functioning  properly. 

The  testing  and  evaluation  of  the  MPDATA  scheme  began  with  the  acquisition  of  the 
FORTRAN-77  code  for  a  one-dimensional  version  of  the  MPDATA  scheme  from  Dr.  Piotr 
Smolaildewicz  at  the  National  Center  for  Atmospheric  Research  (NCAR).  A  driver  test  progr^ 
which  advects  a  square-wave  shaped  mass  perturbation  through  a  domain  ad-infinitum  using 
cyclic  lateral  boundary  conditions  was  written.  The  MPDATA  code  was  then  modified  to 
increase  and  measure  its  computational  efficiency.  The  Smolarkiewicz  code  permits  the  user  to: 
(1)  select  the  numerical  order  of  the  advection  scheme  (from  a  Ist-order  donor  cell  scheme  up  to 
a  3rd-order  MPDATA  scheme);  (2)  specify  whether  a  monotonicity  constraint  should  be 
enforced;  and  (3)  decide  whether  a  specif  correction  for  divergent  flows  should  be  applied.  'The 
performance  of  the  range  of  configurations  was  tested  by  executing  one-dimensional  advection 
experiments  on  an  Apple  Macintosh  II  microcomputer.  The  results  of  the  simple  advection 
experiments  after  5000  timesteps  are  shown  in  Figure  S-3.  It  can  be  seen  that  the  mass  field  was 
severely  diffused  and  the  amplitude  of  the  square  wave  was  reduced  to  46%  of  its  original 
amplitude  when  the  mass  was  transported  with  an  upstream-differencing  (donor  cell)  scheme. 
The  application  of  the  higher  order  MPDATA  scheme  corrected  these  problems.  The  gradients 
of  the  mass  field  were  sharpened  dramatically  by  the  higher  order  formulations.  However,  both 
the  2nd-order  and  3rd-order  versions  of  the  scheme  generated  a  square  wave  with  an  amplitude 
that  was  about  6%  higher  than  its  original  value.  Enforcing  monotonicity  on  the  3rd-order 
solution  eliminated  the  high  frequency  ripples  and  the  overshooting.  The  application  of  this 
constraint  resulted  in  less  than  a  .08%  error  in  the  amplitude  after  50()0  timesteps.  However,  the 
computational  price  for  the  increased  accuracy  was  substantial.  'The  execution  was  2.5  times 
slower  with  the  monotonicity  constraint  than  in  the  case  of  the  3rd-order  scheme  without  the 
monotonicity  constraint.  A  similar  reduction  in  computational  efficiency  was  also  noted  when 
the  accuracy  of  the  numerical  scheme  was  increased  from  2nd-order  to  3rd-order. 
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1-D  MPDATA  Square  Wave  Advection 
(5000  Timesteps) 
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Figure  S-3.  Plot  of  the  mass  field  efter  5000  timesteps  for  the  one-dimensional  square-wave 
advection  aperiments  to  test  the  set  of  configurations  available  in  the  Smolarkiewicz  MPDATA 
scheme. 

After  the  one-dimensional  advecdon  experiments  were  completed,  the  code  for  the  two- 
dimensional  version  of  the  MPDATA  scheme  was  obtained  from  Dr.  Smolarkiewicz  at  NCAR. 
This  code  was  inserted  into  a  new  test  program  which  was  designed  to  advect  a  conical-shaped 
chunk  of  mass  through  an  "infinite"  two-dimensional  grid  domain.  The  two-dimensional 
MPDATA  code  had  the  same  configuration  options  as  the  one-dimensional  code. 

The  initial  conditions  for  the  two-dimensional  advection  experiments  consisted  of  a  cone 
of  mass  with  a  height  of  1  unit  and  a  base  radius  of  5  grid  units  located  near  the  left  boundary  of 
a  30  X  30  cartesian  grid  domain  (Figure  5-4).  The  time  integrations  were  carried  out  to  500 
timesteps  on  an  Apple  Macintosh  Ilex  microcomputer  running  at  a  clock  speed  of  25  MHz.  A 
depiction  of  the  mass  field  at  the  end  of  each  experiment  is  presented  in  Figure  5-5  and  the 
maximum  amplimde  and  relative  computational  efHciency  of  each  scheme  is  listed  in  Table  5-1. 
Figure  5-5a  illustrates  the  severe  degradation  of  the  shape  of  the  mass  penurbation  by  the 
upstream  (donor  cell)  scheme.  In  this  experiment  the  maximum  value  of  the  cone-shaped 
pmurbation  was  reduced  to  37%  of  its  original  value.  The  shape  of  the  cone  collapsed  as  strong 
numerical  diffusion  smeared  the  solution  in  the  direction  of  the  flow.  The  2nd-order  MPDATA 
algorithm  retained  71%  of  the  original  peak  but  it  was  still  somewhat  diffusive  in  the 
alongstream  direction  (Figure  5-5b).  The  best  results  were  obtained  with  the  3rd-order 
MPDATA  scheuie.  Almost  94%  of  the  original  height  remained  at  the  end  of  the  integration. 
(Figure  5-5c).  With  this  version  of  the  scheme  the  gradient  in  the  mass  field  was  quite  sharp  and 
only  a  small  portion  of  the  mass  field  was  diffused  near  the  base  of  the  cone.  The  enforcement  of 
the  monotonicity  constraint  on  the  3rd-order  solution  did  not  improve  the  results.  In  fact,  the 
antidiffusive  flux-limiting  process  "clipped"  the  peak  of  the  cone  (Figure  5-5d).  This  was  not 
totally  unexpected  since  Smolarkiewicz  (1983a)  had  noted  that  the  incorporation  of  a  flux-limiter 
is  extremely  effective  in  preventing  overshooting  when  advecting  shocks  or  contact 
discontinuities  as  in  the  case  of  the  square-wave  used  in  the  one-dimensional  advection 
experiments,  but  it  can  diminish  real  peaks  and  valleys  in  smoothly  varying  distributions. 
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The  "nintiine"  ratio  presented  in  Table  S-1  is  a  measure  of  the  amount  of  CPU  time 
required  for  the  computadons  of  one  of  the  higher  order  versions  of  MPDATA  reladve  to  the  fast 
but  highly  diffusive  donor  cell  scheme.  The  computational  cost  was  increased  by  a  factor  of  4 
and  14  when  the  accuracy  of  the  numerical  scheme  was  increased  to  2nd-order  and  Srd-order 
respectively.  The  execution  was  63%  slower  when  the  inonotonicity  constraint  was  enforced. 
However,  these  comparisons  are  valid  for  scalar  computations  only.  The  results  were  expected 
to  be  significantly  different  for  computers  which  have  the  capability  to  vectorize  the 
computations. 

Table  5-1.  Maximum  amplitude  of  the  mass  field  at  the  end  of  the  two-dimensional  advection 
experunents  and  the  relative  computational  cost  for  each  corfiguration  of  the  MPDATA  scheme 
that  was  tested. 


SCHEME 

MAX  AMPLITUDE 

RUNTIME  RATIO 

Donor  Cell 

.372 

1.00 

2nd-order  MPDATA 

.712 

4.43 

3Td-order  MPDATA 

.939 

14.25 

3rd-order  monotonic  MPDATA 

.864 

23.23 
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Figure  5-4.  Depiction  of  the  initial  conditions  used  in  the  two-dimensional  advection 
experiments  with  the  MPDATA  scheme.  A  conical-shaped  piece  of  mass  is  inserted  near  the 
upstream  boundary  of  the  grid  domain.  Cyclic  lateral  boundary  conditions  are  employed  in  the 
time  integrations. 
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Figure  5-5.  Depiction  of  the  mass  field  after  500  timesteps  with:  (a)  the  upstream  "donor  cell" 
advection  scheme;  (b)  2nd-order  MPDATA  scheme;  (c)  3rd-order  MPDATA  scheme;  and  (d) 
3rd-order  MPDATA  scheme  with  monotonicity  constraint. 
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Before  the  two-dimensional  version  of  the  MPDATA  code  was  installed  in  the  MASS 
model,  two  modifications  were  made  to  the  original  software.  First,  a  new  hybrid  type  of 
boundary  condition  was  in^lemented.  It  consisted  of  the  application  of  a  constant  time-tendency 
for  the  evolution  of  the  mass  fleld  at  the  inflow  boundary  joints  and  an  upstream  or  donor  cell 
advection  scheme  at  outflow  bound^  points.  The  second  modificadon  was  the  incorporation  of 
a  ciffvilinear  coordinate  transformation  "metric”  term,  "m",  into  the  transpon  equadon: 


dy 

dT 


+  m^ 


(5-10) 


to  account  fOT  the  distortions  due  to  the  map  projection  scale  factor  when  the  calculadons  are 
done  on  a  map  image  plane. 

A  series  of  advecdon  experiments  to  test  the  new  boundary  condition  formuladon  was 
executed.  The  inidal  conditions  consisted  of  a  square  block  of  mass  with  a  height  of  1  and  a 
length/width  of  S  grid  units,  located  in  the  center  of  a  30  x  30  grid  domain  (Figure  5-6).  The  3rd- 
order,  monotonic  version  of  MPDATA  was  selected  as  the  transpon  algorithm.  Time  integrations 
were  carried  out  to  300  timesteps  on  an  Apple  Macintosh  Ilex  microcomputer.  The  results  are 
presented  in  Figure  5-7.  The  upstream-differencing  advection  scheme  smoothly  transponed  the 
mass  through  the  outflow  boundary.  Because  this  upstream  scheme  is  much  more  numerically 
diffusive  than  the  antidiffusive  flux-corrected  MPDATA  scheme,  the  leading  gradient  of  the 
mass  field  became  somewhat  degraded.  In  particular,  there  is  a  slight  artificial  "ramping"  in  the 
solution  after  80  timesteps,  when  the  steepest  slope  approaches  the  boundary  (Figure  S-7a). 
Also,  the  leading  edge  of  the  "flat  top"  of  the  block  of  mass  is  prematurely  eroded  at  ^e  outflow 
boundary  gridpoints  (Figures  5-7b  and  5-7c).  Nevertheless,  the  performance  of  this  scheme  was 
deemed  to  be  satisfactop'. 

Before  it  was  inserted  into  the  three-dimensional  MASS  code,  the  two-dimensional 
MPDATA  test  program  was  ported  to  the  Stardent  750  workstation  and  its  computational 
efficiency  was  compared  with  the  performance  on  the  Apple  Macintosh  Ilex  microcomputer. 
The  30  X  30  gridpoint  run  which  required  over  5  minutes  to  execute  on  the  "Mac"  required  only 
12  seconds  on  the  Stardent  750  system  when  the  vectorizaiion  option  was  turned  on. 


Figure  5-6.  Depiction  of  the  initial  conditions  used  in  the  two-dimensional  advection 
experiments  to  test  the  modified  lateral  boundary  condition  and  the  incorporation  of  the 
curvilinear  metric  term  into  the  transport  equation. 
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Figure  S-7.  Depiction  of  the  mass  field  after:  (a)  80  timesteps;  (b)  100  timesteps;  and  (c)  ISO 
timesteps  in  the  simulation  experiment  used  to  test  the  hybrid  inflow-outflow  boundary 
conditions. 
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53  Implementation  of  MPDATA  in  the  MASS  model 

The  originai  MPDATA  code  is  based  upon  the  upstream  advection  scheme  on  a  grid  in 
which  the  advecting  velocities  are  placed  at  points  which  are  located  halfway  between  the  points 

at  which  the  quantity  to  be  advected  (V)  is  deHned.  A  example  of  such  a  grid  is  shown  in  Figure 
S>8.  The  arrangement  of  variables  shown  in  Bgure  S-8  is  commonly  referred  to  as  the  Arakawa 
X"  grid.  However,  the  MASS  model  utilizes  an  unstaggered  (Arakawa  type  "A")  grid  in  which 
all  of  the  prognostic  variables  are  calculated  at  the  same  grid  point  Thus,  the  wind  components 
must  be  interpolated  from  their  unstaggered  positions  to  the  locations  depicted  in  Figure  S-8 
before  the  MPDATA  calculations  can  performed.  Once,  the  interpolation  has  been  completed 
the  remainder  of  the  MPDATA  calculations  can  be  performed  in  the  "C  grid  environment. 
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Figure  S-8.  A  schematic  depiction  of  the  placement  of  variables  on  the  Arakawa  "C"  grid 
utilized  in  the  MPDATA  scheme. 

The  first  step  in  the  MPDATA  calculations  is  to  obtain  an  estimate  of  the  advected 
variable  at  the  next  time  step  from  an  upstream  calculation.  In  the  x-coordinate  direction  the 
calculation  can  be  written  as, 

Vi  =  V"  -  (F  (V?-  V,"  1 '  )-  F  (vf-i  •  VT-  u"i/2  )|  .  (5-11) 

where  the  Fs  are  the  fluxes  of  the  advected  quantity  at  a  location  one  half  grid  increment  to  the 
left  and  right  of  the  point  at  which  the  advection  calculation  is  being  performed.  The  flux  at 
these  points  is  defined  as 


F(Vi.Vi+].u)  =  [(u-i-luDvi  +  (u-|ul)Vi.,j]-^  ,  (5-12) 

2  Ax 

where  At  is  the  time  step  and  Ax  is  the  grid  spacing.  If  were  to  be  used  as  the  updated 
quantity,  the  scheme  would  be  a  pure  upstream  advection  scheme.  Unfonunately,  this  scheme  is 
a  first-order  accurate  scheme  in  both  time  and  space,  and  has  strong  implicit  diffusion  as  shown 
by  Smolarkiewicz  (1983a).  In  order  to  remove  the  implicit  diffusion  in  this  scheme,  MPDATA 
performs  a  second  upstream  advection  step  using  an  anificiaJ  "antidiffusion  velocity"  to  advect 
the  field  in  place  of  the  actual  velocity.  The  "antidiffusion  velocity"  is  defined  as 
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(5-13) 


-  (l«t^i/2l  Ax  -  At  u?  1/2 )  (Vi*i  -  Vi ) 

Ui+i/2  - - ; — ‘-r- 

I  Vi  +  Vi+i  +  e)Ax 

where  e  is  a  small  value  (e.g.  lO*^^)  which  forces  the  antidiffusion  velocity  to  be  zero  when  the 
scalar  quantity  is  zero  at  both  the  i  and  Rl  points.  The  antidiffusion  velocity  is  then  used  to 
advect  the  scalar  field  in  the  second  step  of  the  scheme, 

Vi*  »  V*  -  (f (v*- Vw •  «t+i/2 )•  F (v*-i-  Vi*-  Ui-i/2 )1 .  (5*14) 

and  produce  a  revised  estimate  (V**)  of  advected  quantity  at  grid  point  i.  This  value  can  be 
used  as  the  final  result  of  the  advecdon  process  for  grid  point  i.  If  this  is  the  case  then  this  would 
be  the  2nd-offder  version  of  the  MTOATA  scheme.  However,  the  process  can  be  repeated  again 

by  using  equation  5-13  to  calculate  a  new  antidiffusive  velocity  with  the  v**  values.  This 
procedure  yields  a  higher  order  version  of  MPDATA. 

A  limitCMi  implementation  of  the  MPDATA  scheme  was  made  in  the  MASS  model.  In 
this  implementation,  MPDATA  can  be  utilized  only  for  the  advecdon  of  cloud  water  (ice),  rain 
water  (snow)  or  tracer  material.  In  principle,  the  scheme  could  be  used  for  the  advecdon  of  any 
of  die  prognostic  variables.  However,  in  order  to  assure  dynannical  consistency  among  the 
variables  ^at  are  strongly  coupled,  the  scheme  must  be  used  with  all  of  these  variables 
simultaneously.  If  it  is  implemented  with  only  one  variable  (for  example,  the  water  vapor 
mixing  ratio)  the  slight  differences  in  phase  spi^  between  the  MPDATA  scheme  and  oAer 
schemes  such  as  Adams-Bashforth  can  result  in  the  generation  of  dynamically  inconsistent  fields 
and  the  generation  of  spurious  features. 

Once  the  installation  of  the  MPDATA  subroutine  into  the  MASS  model  was  completed,  a 
series  of  simple  simulations  were  executed  to  verify  that  the  code  was  correctly  installed  and  to 
test  the  performance  of  the  scheme.  The  first  experiment  was  the  advecdon  of  a  massless  tracer 
with  an  idealized  wind  distribution  on  a  SO  km  grid.  In  this  test,  the  wind  was  constrained  to  be 
25  m  s*^  from  the  west  at  all  model  grid  points  throughout  the  simulation.  That  is,  the  wind  was 
spatially  uniform  and  time  invariant.  The  left  column  illustrates  the  total  mass  in  kg  m*^  that  was 
instantaneously  injected  into  the  model  atmosphere  at  the  initial  time  (1200  UTC).  This  value 

was  the  result  of  instantaneously  injecting  1.0  x  10^^  kg  of  tracer  mass  per  model  grid  cell  per 
500  meter  vertical  layer.  As  the  illustrations  in  the  left  column  of  Figure  5-9  suggest,  the  tracer 
substance  was  injected  into  only  4  MASS  model  colutims  over  southern  California.  The 
underlying  geography  in  this  depiction  is  used  for  reference  purposes  only.  There  was  no  terrain 
or  other  geography-related  effects  incorporated  into  the  simulation.  The  use  of  only  4  horizontal 
points  poses  a  severe  test  for  an  advecdon  scheme  since  the  lateral  scale  of  the  tracer  pool  is  at 
the  lower  limit  of  what  is  resolvable  on  the  grid  system.  The  top  row  of  Figure  5-9  (a-c) 
illustrates  results  from  the  2nd-order  version  of  MPDATA,  the  middle  row  (d-f)  depicts  the 
performance  of  the  3rd-order  version  of  MPDATA  while  the  bottom  row  (g-i)  shows  the 
performance  of  the  4th-order  accurate  space  difference/Adams-Bashforth  time  integration 
advecdon  scheme  (the  miginal  MASS  advecdon  scheme).  In  this  case,  a  good  performance  by  an 
advecdon  scheme  would  result  in  the  translation  of  the  injected  mass  downwind  while  preserving 
its  shape  and  maximum  concentration  value.  Figure  5-9  illustrates  that  all  three  of  the  schemes 
reduce  the  peak  value  significantly  during  the  6  hour  simulation  period.  Surprisingly,  the  original 
MASS  scheme  produces  a  concentration  value  that  is  closest  to  Ae  value  at  the  time  of  the  initial 
injection,  llie  MPDATA  scheme  had  been  expected  to  preserve  the  shape  and  maximum  value 
of  the  mass  pool  better  than  the  original  MASS  scheme.  These  surprising  results  are  most  likely 
dependent  on  the  scale  of  the  tracer  pool.  Apart  from  this  somewhat  surprising  outcome,  the 
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results  from  these  test  simulations  confirmed  the  expectations  that:  (1)  the  3rd-oider  MPDATA  is 
better  at  preserving  the  maximum  value  of  the  tracer  pool  than  the  2nd-order  scheme.  (2)  the 
original  MASS  scheme  is  much  noisier  than  the  MPDATA  scheme;  and  (3)  as  advertU^  the 
positive-definite  MPDATA  scheme  produced  no  negative  values  of  tracer  mixing  ration  while 
the  miginri  MASS  scheme  produced  some  negative  values  in  the  vicinity  of  large  tracer  mixing 
ratio  g^ents. 

The  damping  of  the  tracer  pool  in  these  initial  experiments  was  examined  more  closely.  It 
was  (tetermined  that  the  scheme  worked  well  when  running  advecdon  experiments  with  an  initial 
tracer  field  consisting  of  a  square/sine  wave  resolved  by  4  or  more  grk^mints  in  the  direction  of 
the  flow.  However,  the  inidal  distribution  of  tracer  used  in  the  initial  experiments  conducted  with 
the  MASS  model  was  only  1  to  2  gridpoints  wide  in  the  direction  of  the  flow. 

This  behavitx'  was  traced  to  the  fact  that  numerical  finite  differencing  schemes  which  use 
flux  correction  techniques  (MPDATA  and  FCT,  for  example)  calculate  their  flux  limiters  based 
on  the  minimum  and  maximum  values  at  neighboring  midpoints.  That  is,  they  decide  how  much 
additional  mass  should  be  added  to  the  diffused  solution  at  a  gridpoint  in  order  to  restore  the 
mass  field  to  its  pre-diffused  form  by  looking  at  peaks  and  valleys  at  surrounding  gridpoints.  If  a 
“spike”  of  tracer  is  not  adequately  resolved,  then  the  extrema  can  “fall  through  the  grid  mesh” 
during  the  advecdon  process,  never  to  be  recovered  by  the  flux  limiters  again.  Hiis  numerical 
“clipping"  process  can  often  be  remedied  by  using  a  “smarter”  flux-limidng  algorithm,  one 
which  “remembers"  the  values  of  the  extrema  for  the  previous  timestep  or  one  which  calculates 
the  minimum  and  maximum  values  at  higher  than  grid  mesh  resoludon.  This  type  of  approach 
was  invesdgated  by  coding  a  version  of  Zalesak’s  (1979)  flux  limiter,  which  calculates  the 
maxima  and  minima  between  grid  points,  by  finding  the  intersection  of  two  lines  which 
approximate  the  shape  of  the  mass  field  at  half-gridpoints.  A  series  of  one-dimensional  advecdon 
experiments  was  then  made  to  assess  the  performance  of  the  new  flux  limiter.  The  results  of 
these  experiments  are  shown  in  Figure  5-10.  A  one  gridpoint  wide  disturbance  of  amplitude  1.0 
was  insened  into  the  inidal  mass  field  and  transported  with  a  constant,  uniform  velocity.  The 
differtiices  are  quite  apparent  even  after  only  25  dmesteps.  The  new  Zalesak  flux  limiter, 
implemented  in  the  flux-corrected  transpon  (FCTT)  algorithm,  recovered  about  37%  of  the 
original  amplitude,  while  the  old  limiter  only  retained  31%.  The  width  of  the  mass  field  is  also 
narrowed,  which  results  in  a  more  realistic,  “spike-like”  representadon.  The  amplitude  sdll  falls 
short  of  the  Adams-Bashforth  scheme,  which  holds  onto  59%  of  the  original  peak  value. 
However,  it  should  be  noted  that  the  two  FCT  soludons  are  much  less  noisy  in  the  originally 
unperturbed  regions  than  the  Adams-Bashforth  solution.  The  latter  advecdon  algorithm  leaves  a 
long  wake  of  trailing  ripples  (numerical  noise  or  Gibbs  oscilladons)  due  to  numerical  wave 
dispersion,  with  amplitudes  as  high  as  37%  of  the  original  height. 

Dr.  Piotr  Smolarldewicz  was  consulted  k>  see  if  he  h^  developed  a  better  algorithm  for 
transporting  such  “spike-like”  distribudons  of  tracer.  He  indicated  that  he  had  tried  customized 
flux  limiters  in  the  past,  but  had  concluded  that  once  there  is  a  very  narrow  spike  in  a  field,  finite 
difference  methods  are  simply  inadequate.  He  indicated  that  the  reason  for  this  can  be  understood 
by  expanding  the  finite  difference  form  of  the  governing  advection  equation  in  the  form  of  a 
Taylor  series.  He  indicated  that  the  result  be  that  the  expansion  will  approximate  the  original 
analydcal  equadon  plus  some  extra  terms  (i.e.  error  terms).  The  error  terms  scale  as  NP,  where  N 
is  the  number  of  grid  points  per  scale  of  interest  and  P  is  the  order  of  the  term.  Therefore,  no 
matter  how  accurate  the  Kheme  is.  if  the  number  of  ^d  points  per  scale  of  interest  is  1,  the 
Reynolds-like  number  will  be  1  too.  In  other  words,  if  a  one-gridpoint  perturbadon  is  to  be 
simulated,  even  a  super  accurate,  say,  lOOth-order  scheme  will  suffer  tom  errors  comparable  to 
the  analydc  terms.  Thus,  finite-difference  techniques  are  not  appropriate  for  under-resolved 
features.  There  are  techniques  which  claim  the  ability  to  treat  under-resolved  features,  for 
example  moments  methods.  But  all  they  do  is  store  more  variables  (values  per  gridpoint)  which 
is,  in  essence,  an  increase  in  the  resoludon  of  the  grid. 
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Furthermore,  the  Zalesak  flux  limiter  is  customized  to  handle  isolated  extrema.  It  creates 
artiflcial  ovenhooting  when  the  initial  tracer  is  more  unifonn  in  shape,  as  shown  in  Figure  5-11. 
In  general,  specialized  flux  limiters  should  not  be  implemented  unless  the  shape  of  the  initial 
tracer  field  is  known  beforehand.  This  is  clearly  not  the  case  in  the  general  application  of  a 
three-dimensional  mesoscale  model. 

Advection  experiments  were  conducted  with  more  realistically  shaped  tracer  distributions 
in  order  to  evaluate  the  relative  accuracy  of  the  FCT  and  MPDATA  schemes.  Transpon 
simulations  using  initial  square  wave  and  sine  wave  distributions  all  produced  similar  results 
with  the  MPDATA  scheme  retaining  more  of  the  initial  disturbance  amplitude  and  preserving 
steep  gradients  better  than  flux-cmected  transpon  (Figure  5-12).  The  higher  numerical  diffusion 
in  the  FCT  algorithm  became  especially  apparent  in  runs  which  extended  beyond  1(XX)  timesteps. 

These  results  suggested  that  the  ^DATA  scheme  without  the  Zalesak  flux  limiter  was 
the  best  choice  for  use  in  the  MASS  model.  The  scheme  was  then  configured  to  advect  the  water 
vapor  mixing  ratio  in  the  three-dimensional  version  of  the  MASS  m^l.  In  order  to  test  the 

scheme  in  a  complex  wind  pattern  the  model  was  initialized  with  a  data  matrix  of  85  x  55  x  15 
points  from  1200  UTC  AprU  10, 1979  (the  day  of  the  infamous  Wichita  Falls  tornado).  Three  6- 
hr  simulations  were  executed  from  this  initial  dataset.  In  order  to  have  a  pure  comparison  of  the 
performance  of  each  of  the  horizontal  advection  schemes,  all  other  processes  (e.g.  vertical 
advection  of  vapor,  condensation,  etc.)  which  act  on  the  water  vapor  field  were  omitted  from  the 
simulation.  The  first  simulation  employed  the  4th-oider  accurate  centered  space  difference  and 
the  Adams-Bashforth  time  integration  scheme  (i.e.  the  original  MASS  advection  scheme)  to 
horizontally  advect  water  vapor.  The  second  simulation  was  identical  to  the  first  except  that  the 
donor  cell  version  of  MPDATA  was  used  for  the  horizontal  moisture  advection.  The  third 
simulation  used  the  3rd-order  version  of  MPDATA  with  the  monotonicity  constraint  and  a 
correction  for  divergent  flow.  The  dew  point  for  the  lowest  model  layer  (approximately  30  mb 
above  the  surface)  for  each  of  the  simulations  is  shown  in  Figure  5-12.  A  comparison  of  the  field 
produced  by  the  either  of  the  MPDATA  schemes  (Figures  5- 12b  and  5- 12c)  with  that  produced 
by  the  Adams-Bashforth  scheme  (Figure  5- 1 2a)  immediately  reveals  that  MPDATA  produces 
much  sharper  gradients  and  maxima  and  minima  that  are  more  extreme.  In  fact,  in  certain  areas 
the  extreme  values  appear  to  be  slightly  excessive.  This  is  probably  due  to  the  absence  of  other 
processes  which  offset  the  advection  effects  in  a  full  physics  simulation.  Another  difference 
between  the  simulations  can  be  seen  by  examining  the  dew  point  pattern  over  southwestern 
Texas  and  adjacent  Mexico.  In  this  area,  the  Adams-Bashforth  simulation  (Figure  5- 12a) 
produces  noticeable  high-frequency  oscillations  in  proximity  to  the  large  moisture  gradient  over 
Texas  while  no  trace  of  oscillations  can  be  found  in  either  of  the  MPDATA  simulations.  A 
second  simulation  which  utilized  the  complete  moisture  physics  was  also  executed.  In  this 
simulation  the  temperature  was  advected  with  the  Adams-Bashforth  scheme  and  the  water  vapor 
was  advected  with  the  3rd  order  MPDATA  scheme.  This  simulation  generated  spurious  areas  of 
saturation,  latent  heating  and  precipitation  due  to  the  differences  in  numerical  dispersion  and 
phase  speeds  between  the  two  advection  schemes.  This  experiment  supported  the  notion  that  it  is 
not  wise  to  utilize  different  advection  schemes  for  prognostic  variables  that  are  strongly  coupled. 

In  order  to  test  the  ability  of  the  MPDATA  scheme  to  advect  a  mass  of  particles  injected 
into  the  atmosphere  in  a  real-world  situation,  a  three-dimensional  simulation  of  the  ash  cloud 
generated  by  the  eruption  of  Mt.  St.  Helens  on  May  18,  1980  was  executed.  A  version  of  MASS 

which  used  a  90  x  60  x  17  matrix  and  a  15  km  grid  increment  was  used  for  this  experiment.  The 
vertical  profile  and  rate  of  injection  of  ash  was  estimated  from  data  published  by  Sama-Wojcicki 
et  al.  (1982)  and  Carey  and  Sigurdsson  (1982).  A  3rd-order  version  of  MPDATA  without  any 
lateral  numerical  diffusion  was  used  for  this  simulation  experiment.  A  depiction  of  the  time 
evolution  of  the  vertically  integrated  (surface  to  the  top  of  the  model  domain)  airborne  ash  for  the 
8  hr  period  after  the  initial  eruption  is  shown  in  Figure  5-13.  The  bold  dashed  lines  in  this  figure 
indicate  the  outline  of  the  satellite-observed  ash  cloud  at  each  time  as  diagnosed  by  Sama- 
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Wojtckt  et  al.  (1982).  The  transpon  simulated  by  the  MASS  model  with  the  MPDATA 
advection  scheme  seems  quite  reasonable.  The  small  discrepancies  between  the  simulated  ash 
distribution  and  the  satellite-observed  ash  cloud  pattern  are  probably  attributable  to:  (1)  errors  in 
the  model  wind  forecast  because  of:  (a)  the  customary  limitadons  in  the  amount  arui  distribudon 
of  inidalizadon  data,  and  (b)  limitations  in  model  resolution  and  the  physics;  (2)  imprecise 
specificadon  of  the  verdcal  prorile  and  rate  of  ash  injection  into  the  amrosphere  from  the 
erupdtm;  and  (3)  the  unknown  value  of  the  vertically  integrated  ash  content  that  corresponds  to 
the  boundary  of  the  ash  cloud  in  die  satellite  imagery. 
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Figure  5-9.  The  concentration  (kg  nv^)  of  a  massless  tracer  integrated  from  the  surface  to  609  m 
(2(XX)  ft)  from  three  advection  experiments.  The  winds  are  spatially  uniform  and  time  invariant 
at  a  value  of  25  m  s  f  Tracer  distributions  at  the  initial  time  (1200  UTC).  3  hours  (1500  UTC) 
and  6  hours  (1800  UTC)  after  initialization  are  shown  for  each  experiment,  (a-c)  2nd-order 
MPDATA;  (df)  3rd-order  MPDATA;  (g-i)  Adams-Bashforth  with  fourth  order  spatial  finite 
differencing. 
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Spike  Advcction  Experiment:  25  Timesteps 


Figure  S-IO.  The  results  from  a  series  of  one-dimensional  advection  experiments  executed  to  test 
the  performance  of  the  Zalesak  one-point  extremum  flux  limiter.  A  4th-order  version  of  the  flux- 
corrected  of  the  flux-corrected  transport  (FCT)  algorithm  was  used  to  transport  a  one-gridpoint 
"spike"  of  tracer  with  an  amplitude  of  1.0  throu^  a  uniform  flow  field.  The  simulations  show 
that  the  Zalesak  limiter,  denoted  by  NL  for  "new  limiter",  is  able  to  retain  more  of  the  initial 
arriplitude  of  the  disturbance  and  replicate  the  narrowness  of  the  spike-like  representation  more 
faitlrfully  than  the  old  flux  limiter  (OL).  The  results  of  he  4th-order  Adams-Bashforth  scheme, 
which  is  currently  used  in  the  MASS  model,  are  also  presented.  Cyclic  boundary  conditions 
were  used  for  all  experiments. 

Square  Wave  Advection  Experiment: 

300  Timesteps 
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Figure  5-11.  The  results  from  of  a  1-D  square -wave  advection  experiment  which  show  how  the 
Zalesak  one-point  extremum  flux  limiter  (NL  for  "new  limiter"),  designed  to  recover  unresolved 
minima/maxima,  "overshoots"  the  amplitude  of  the  initial  tracer  field  when  it  is  uniform  in 
shape.  Also  notice  the  long  trailing  wake  of  Gibbs  oscillations  (noise)  due  to  numerical  wave 
dispersion  present  in  the  Adams-Baslforth  scheme,  which  is  currently  used  in  the  MASS  model. 
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Figure  5-12.  Dew  point  in  the  lowest  model  layer  from  three  six-hour  MASS  simulations, 
initialized  at  1200  UTC  April  10, 1979.  The  contour  interval  is  4  U.  All  of  the  sirrmlations  are 
identical  with  the  exception  that  the  advection  of  water  vapor  mixing  ratio  is  computed  by  :  (a) 
the  fourth  order  space  difference  /  Adams-Bashforth  time  difference  scheme,  (b)  the  donor  cell 
version  of  MPDATA,  and  (c)  the  3rd-order  MPDATA  scheme  with  monotonicity  constraint  and 
correction  for  divergent  flow. 


SBIR  Phase  II  Final  Report 


DAAD07-90-C-0134 


Page  46 


Figure  S-J3.  Vertically-integrated  airborne  ash  per  unit  area  (kg  m-^)  from  a  15  km  MASS 
simulation  with  the  3rd-order  MPDATA  scheme  used  for  ash  transport.  The  ash  depictions  are 
for:  (a)  1745  UTC;  (b)  1945  UTC:  (c)  2145  UTC:  and  (d)  2345  UTC.  The  first  contour  at  each 
time  is  .05  kg  m'^  ar^  the  contour  interval  is  05  kg  m'^.  The  bold  dashed  lines  indicate  the 
outline  of  the  satellite-observed  ash  cloud  at  each  time  as  reported  by  Sarna-Wojicki  et  al. 
(1982). 


SBIR  Phase  n  Final  Report 


DAAD07-90-C-0134 


Page  47 


6.  Data  Assimilation  System 

The  initimlization  of  mesoscale  numerical  models  in  a  battlefield  environment  is  a 
significant  challenge.  In  such  an  environment,  observatitmal  data  is  typically  available  from  a 
variety  of  sources,  each  with  its  own  strengths  and  weaknesses.  Any  single  data  source  may  be 
available  in  mily  a  portion  of  the  area  of  interest,  may  not  measure  all  ^pendent  atmospheric 
vaiiaUes  or  may  contain  individual  observations,  each  valid  at  a  different  time.  The  challenge  is 
to  combine  the  information  availdt>le  from  each  data  source  into  a  well-balanced  m^el 
initialization  that  includes  all  mesoscale  features  which  are  resolvable  by  the  model. 

In  mder  to  take  the  best  possible  advantage  of  the  various  types  of  data,  modelers  have 
developed  various  methods  of  continuously  inserting  current  and  past  data  into  a  numerical 
model  during  a  q)e^ed  pre-fmecast  period.  This  method,  now  known  as  four  dimensional  data 
assiimladon  (FDDA),  allows  the  model  inidalizatimi  to  take  maximum  advantage  of  data  sources 
available  at  asynopde  times  as  weU  as  data  with  hi^  temporal  resolution.  FDDA  was  pioneered 
by  Chamey  et  al.  (1969)  who  suggested  that  the  model's  prognostic  equations  would  provide 
time  continuity  and  dynamic  coupling  between  the  various  fields. 

The  MASS  model  uses  several  types  of  FDDA  and  an  innovative  method  of  static 
moisture  initialization,  in  addition  to  the  standard  static  objective  analysis  with  a  fust  guess. 
These  include:  (1)  Newtonian  relaxation  or  nudging  of  gridded  rawinsonde,  surface  or  profiler 
obsovadons  (Hoke  and  Anthes,  1976)  in  which  &e  model  state  is  relaxed  towards  an  analysis  of 
observations;  (2)  the  use  of  Manually  Digitized  Radar  (MDR)  data  to  specify  moisture 
convergence  (precipitation  rates)  in  areas  which  are  (are  not)  subjea  to  convection  according  to 
the  Kuo-MESO  cumulus  parameterization  scheme;  and  (3)  the  enhancement  of  the  three 
dimensional  moisture  analysis  with  surface  observations  of  clouds,  pilot  reports,  MDR  data  and 
infrared  satellite  images. 

6.1  Surface  Data 

The  assimilation  of  surface  data  is  based  on  the  techniques  of  Stauffer  et  al.,  (1991).  The 
scheme  is  described  in  depth  in  the  MASS  Reference  Manual  (MESO,  1993a).  It  was  initially 
implemented  in  a  manner  identical  to  that  of  Stauffer  et  al.,  (1991).  As  the  scheme  was  tested 
though,  several  improvements  were  made.  These  include: 

1)  The  use  of  an  Ekman  wind  profile  to  nudge  winds  in  the  planetary  boundary 
layer  (PEL)  in  regions  where  a  profile  can  be  found  that  achieves  a  reasonable 
fit  with  both  the  gridded  surface  wind  analysis  and  the  model  winds  at  the  top 
of  the  PEL.  This  is  superior  to  assuming  a  constant  wind  velocity  in  the  PEL. 

2)  The  adjustment  of  gridded  surface  wind  analysis  to  reflect  the  model-defined 
surface  roughness  at  each  individual  grid  point.  This  was  done  mainly  in 
response  to  the  tendency  of  surface  nudging  to  spuriously  reduce  wind  speeds 
over  water,  since  most  observations  are  land-based. 

3)  The  nudging  of  gridded  surface  temperature  in  the  PEL.  Teiiqierature  is  nudged 
most  strongly  in  the  lowest  model  layer.  The  strength  of  nudging  decreases 
until  it  reaches  0  at  the  top  of  the  PEL.  Stauffer  et  al.,  (1991)  did  not  nudge 
surface  temperature  because  of  its  potentially  harmful  interaction  with  the  PEL 
parameterization  (i.e.,  the  PEL  parameterization  will  be  forced  into  another 
regime).  Surface  temperature  nudging  in  the  MASS  model,  however,  has 
pi^uced  mostly  positive  effects.  Perhaps  this  is  due  to  the  fact  that  errors  in 
the  representation  of  cloudiness  along  with  incomplete  knowledge  of  the 
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distribution  of  soil  types,  soil  moistuze  and  vegetation  can  sometimes  lead  to 
significant  temperature  errors  in  the  model  PBL.  While  nudging  temperature 
will  not  prevent  these  limitations  from  leading  to  temperature  errors  after  the 
forecast  period  has  begun,  it  can  minimize  the  propagation  of  PBL  temperature 
errors.  Another  reason  for  nudging  temperature  is  that  in  regions  with  large 
errors  in  the  simulated  temperature  Held,  the  model  PBL  is  likely  to  be  in  a 
different  regime  than  the  actual  atmosphere  anyway. 

4)  Numerous  modifications  were  made  to  the  Barnes  objective  analysis  scheme  to 
improve  the  analysis  in  the  vicinity  of  coastlines.  Coastlines  are  frequently 
collocated  with  sharp  discontinuities  in  low-level  wind,  temperature  and 
moisture  fields.  Most  objective  analysis  schemes  do  not  resolve  these 
discontinuities.  This  problem  is  especially  compounded  by  the  fact  that  the  vast 
majority  of  surface  observations  are  land-based.  With  a  simple  objective 
analysis  scheme,  the  values  of  surface  variables  are  often  heavily  influenced  by 
land-based  observations  at  grid  points  which  are  located  just  offshore.  The 
improved  Barnes  scheme  eliminates  this  problem.  See  Section  7.1  of  the 
f^SS  Reference  Manual  (MESO,  1993a)  for  details  on  the  modifications 
made  to  the  Bames  scheme. 

5)  The  method  of  calculating  confidence  factors  was  improved  so  that  the 
confidence  factor  would  fall  sharply  to  0  where  there  were  discontinuous 
changes  in  data  density  such  as  at  a  coastline.  The  differentiation  between  land 
and  water  observations  also  contributes  to  sharp  changes  in  the  confidence 
factor  at  coastlines.  See  Section  7.1  of  the  MASS  Reference  Manual  (MESO, 
1993a)  for  details. 


6.2  Rawinsonde  Data 

The  assimilation  of  rawinsonde  data  is  described  in  depth  in  the  MASS  Reference 
Manual  (MESO,  1993a).  As  for  the  surface  data,  the  rawinsonde  nudging  scheme  was  initially 
implemented  in  a  manner  identical  to  that  of  Stauffer  et  al.,  (1991).  As  the  scheme  was  tested 
though,  several  improvements  were  made.  These  include: 

1)  Changes  in  the  method  of  calculating  the  confidence  factor  so  that  it  would  fail 
sharply  to  0  on  the  sparse-data  side  of  discontinuities  in  data  density.  The 
confidence  factor  is  also  calculated  at  each  level  in  the  vertical  and  for  each 
variable.  Before  these  changes  were  implemented,  the  nudging  of  ^dded  data 
just  on  the  sparse-data  side  of  a  data  discontinuity,  or  in  other  regions  of  low 
data  density,  often  caused  large  oscillations  in  the  mass  field. 

2)  The  nudging  of  rawinsonde  data  only  within  4  hours  of  a  rawinsonde 
observation  time  and  nudging  with  full  weight  only  within  2  hours  of  an 
observation  time.  Rawinsonde  data  is  available  only  every  12  hours  which  is 
too  infrequent  to  resolve  the  finer  scale  structure  that  a  mesoscale  model  can 
simulate.  For  this  reason,  a  linear  interpolation  in  time  is  not  sufficient  mid¬ 
way  between  observation  times. 


6.3  Profiler  Data 

The  assimilation  of  profiler  data  is  based  on  the  techniques  used  by  Stauffer  et  al.,  (1991) 
to  nudge  rawinsonde  data.  Height-based  wind  profiler  observations  are  interpolated  to  pressure 
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surfaces  by  using  a  recent  model  initialization  file  to  estimate  the  height  field.  Temperatures  are 
derived  through  the  inversion  of  the  divergence  equation  as  in  Cram  et  al.,  (1991).  The  scheme  is 
described  in  depth  in  the  MASS  Reference  Manusd  (MESO,  1993a). 

Wind  i^filer  data  provides  accurate  and  fluent  3*dimensional  data  that  seems  to  be 
ideal  for  assimilation  into  a  mesoscale  model.  However,  there  are  several  drawbacks  including  a 
lack  of  direct  infcmnation  on  the  mass  field,  a  lack  of  any  data  within  about  1  km  of  the  earth's 
surface,  and  poor  spatial  resolution  relative  to  the  temporal  resolution.  The  development  of  the 
profiler  nudging  scheme  focused  on  efforts  to  overcome  these  drawbacks. 

The  first  challenge  was  to  interpolate  the  wind  information  to  pressure  coordinates  so  that 
it  could  be  analyzed.  This  was  done  by  assuming  that  a  MASS  output  file  at  or  near  the  profiler 
observation  time  would  provide  a  reasonably  accurate  height-pressure  relationship.  Once  this 
was  accomplished,  the  next  step  was  to  derive  some  sort  of  mass  field  information  from  the  wind 
field,  since  accurate  mesoscale  forecasts  typically  require  an  accurate  representation  of  both  the 
mass  and  the  wind  field.  Cram,  et  al..  (1^1)  invert^  the  divergence  ^nation  to  derive  the  3- 
dimensional  temperature  Held  when  only  the  wind  field  was  known.  This  technique  was  adapted 
to  the  MASS  m^el,  but  a  few  changes  needed  to  be  made  to  the  nudging  scheme.  First,  MASS 
output  was  used  to  provide  a  first  guess  for  the  profiler  wind  analysis.  This  produced  a 
reasonable  wind  field  throughout  the  ni^el  domain  and  prevented  the  divergence  technique  from 
producing  ridiculous  temperatures  in  and  near  data  sparse  regions.  The  MASS  output  was  also 
used  to  provide  the  height  field  at  the  domain  boundarfes.  which  the  technique  requires.  Second, 
the  time  derivative  term  was  ignored  (term  H  in  C!ram  et  al.'s  equation  (1)).  This  term  had  little 
effect  on  the  temperature  derivation,  and  ignoring  it  simplified  the  process  by  eliminating  the 
need  for  a  second  profiler  wind  analysis.  Finally,  tiie  temperature  confidence  factor  was  set  to  0 
within  1750  m  of  the  earth's  surface.  This  was  done  because  it  was  discovered  that  frictional 
effects  which  are  not  included  in  the  divergence  equation,  but  are  present  in  the  PBL, 
contaminate  the  temperature  analysis  at  lower  levels. 

6.4  Manually  Digitized  Radar  (MDR)  Data 

MDR  data  and  reports  of  areal  coverage  of  radar  echoes  are  used  to  estimate  the 
precipitation  rate  at  each  model  grid  point.  This  data  is  then  inserted  into  the  model  by  specifying 
a  parabolic  latent  heating  profile  at  grid  points  which  are  ineligible  for  convection,  according  to 
the  Kuo-MESO  cumulus  parameterization  scheme  or  by  specifying  the  moisture  supply 
parameter  at  points  which  are  eligible  for  convection.  The  scheme  is  described  in  depth  in  the 
MASS  Reference  Manual  (MESO,  1993a).  The  implementation  of  the  MDR  data  assimilation 
scheme  is  described  below. 

The  development  of  the  MDR  assimilation  scheme  went  through  several  phases.  First, 
VIP  levels  and  reports  of  areal  coverage  were  used  to  estimate  cloud  top  levels  and  mean  relative 
humidities.  The  RH  information  was  then  nudged  into  the  model  simulation.  Hiis  scheme 
performed  reasonably  well  in  regions  where  stable  grid  scale  precipitation  was  falling,  but  not  so 
well  where  convection  was  occurring.  The  next  phase  was  to  derive  precipitation  rate  estimates 
from  the  MDR  and  areal  coverage  data.  Estimates  were  fine  tuned  through  trial  and  error.  Then, 
rather  than  nudging  RH  in  convective  regions,  the  precipitation  rate  estimate  was  used  to  specify 
the  moisture  convergence  field  in  the  Kuo-MESO  convective  parameterization  scheme.  This 
produced  a  significant  improvement  in  the  representation  of  convection  during  the  data 
assimilation  period.  The  final  phase  was  in  response  to  the  fact  that  nudging  RH  in  regions  of 
precipitation  was  ineffective  in  regions  where  the  model  incorrectly  forecasted  downward 
vertical  motions.  This  often  occurred  where  mesoscale  precipitation  regions  were  not  resolved  by 
the  observing  network.  To  help  improve  the  simulation  in  such  regions,  the  precipitation  rate  and 
cloud  top  height  estimates  were  used  to  specify  a  parabolic  latent  heating  profile  in  regions  of 
grid  scale  precipitation.  This  was  done  in  place  of  nudging  RH.  Below  an  assumed  cloud  base  of 
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about  1000  m  above  the  ground,  RH  was  nudged  as  it  had  been.  This  final  change  resulted  in  a 
substantial  improvement  in  the  representation  of  grid  scale  precipitation  areas  that  were  not  well 
resolved  by  the  observing  networic. 

Experience  has  shown  that  the  assimilation  of  MDR  data  helps  to  accurately  define  both 
the  grid  si^e  and  subgrid  scale  precipitation  field  at  the  time  of  nuxlel  initialization.  The  benefits 
of  MDR  assimilation  typically  last  for  about  6  to  9  hours.  They  are  most  dramatic  during  the 
initial  three  hours  due  to  the  elimination  of  the  precipitation  spin-up  problem  that  is  common  to 
all  numerical  models. 


6.5  Synthetic  Relative  Humidity 

The  synthetic  RH  retrieval  scheme  uses  visual  observations  of  clouds,  infrared  (IR) 
satellite  data  and  MDR  data  to  enhance  the  rawinsonde  moisture  analysis.  The  scheme  is 
described  in  depth  in  the  MASS  Reference  Manual  (MESO,  1993a).  The  implementation  of  the 
synthetic  RH  analysis  scheme  is  described  below. 

The  main  challenge  in  tiie  implementation  of  the  synthetic  RH  scheme  was  to  determine 
just  what  information  about  tiie  RH  field  could  be  gleaned  from  the  various  data  sources  and  how 
the  different  data  sources  could  be  combined  to  provide  the  best  estimate  of  RH. 

Tlie  derivation  of  vertical  RH  profiles  from  visual  observations  of  clouds  underwent 
several  changes  during  development.  First,  the  mean  RH  at  each  level  above  the  ground  was 
simply  averaged  from  all  of  the  soundings  associated  with  each  observation  in  a  given  cloud 
weather  category.  Next,  statistical  correlations  were  determined  based  on  the  height  and  extent  of 
coverage  of  each  cloud  layer.  Finally,  statistical  correlations  were  stratified  into  layers  so  that  a 
different  set  of  regression  equations  existed  for  levels  that  were  below  the  lowest  cloud  layer, 
between  the  lowest  and  the  next  highest  cloud  layer,  between  the  middle  and  the  highest  layer, 
and  above  the  highest  cloud  layer.  Each  of  these  improvements  helped  to  reduce  the  root  mean 
square  errors  of  the  RH  estimates.  This  resulted  in  more  statistical  soundings  with  a  small 
enough  error  to  be  included  in  the  RH  analysis.  Whenever  relative  humidity  fields  derived  from 
visusd  cloud  observations  are  available,  they  are  blended  with  the  model  RH  analysis  through  an 
objective  analysis  before  any  other  adjustments  are  made.  MDR  and  satellite  data  are  then  used 
to  provide  the  finer  structure  of  the  RH  field. 

Several  data  sources  were  combined  to  determine  the  best  possible  distribution  of  cloud 
coverage,  base  and  height.  IR  satellite  data  is  ideal  for  determining  cloud  top  height  and  areal 
coverage  of  clouds.  Cloud  base  information  can  be  determined  by  visual  observations.  The 
presence  of  MDR  echoes  or  the  detection  of  a  convective  tower  by  Adler  and  Negri's  (1988)  IR 
satellite  technique  indicates  that  a  moist  layer  likely  extends  from  cloud  top  to  the  ground.  When 
satellite  data  is  missing,  MDR  data  can  also  be  used  to  estimate  cloud  top  height  and  areal 
coverage,  as  is  done  in  the  MDR  data  assimilation  scheme  (see  Section  14.5  in  the  MASS 
Reference  Manual,  MESO,  1993a).  Since  IR  satellite  data  has  higher  resolution  and  provides  a 
more  reliable  estimate  of  cloud  coverage,  it  overrides  MDR  data  in  the  determination  of  cloud 
top  height  and  coverage.  When  an  MDR  echo  is  not  present,  it  is  impossible  to  determine 
whether  or  not  a  cloud  layer  is  continuous  between  the  satellite  determined  cloud  top  and  the 
visually  determined  cloud  base.  Experimentation  shows  that  when  both  the  satellite  and  the 
visual  observation  indicate  extensive  cloud  coverage,  it  is  wise  to  err  on  the  side  of  over¬ 
moistening  the  atmosphere  and  assume  the  cloud  is  continuous  in  the  vertical.  When  either  a 
cloud  base  or  cloud  top  height,  but  not  both  can  be  determined,  the  cloud  is  assumed  to  be  thin. 
After  extensive  trial  and  error,  the  synthetic  RH  scheme  evolved  into  the  quite  complicated,  but 
effective  algorithm  that  is  described  in  Chapter  5  of  the  MASS  Reference  Manual  (MESO, 
1993a). 

The  synthetic  RH  scheme  is  most  valuable  when  simulating  cases  with  weak  or 
nonexistent  synoptic  forcing  and  strong  boundary  layer  forcing.  In  these  situations,  convection  is 
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nftMi  triflMred  bv  diffeteatUl  heatiag  due  to  ineso*p  scele  viiiMiwii  in  cloud  C(^».  The 
8ynihelic*fS?8ctone  is  the  component  of  the  model  dau  asshniUiion  system  which  is  most 
Active  at  adding  such  fine  scale  variations  in  cloudiness. 
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7 .  The  Observing  System  Simulation  Experiment 
(OSSE) 

The  data  assimilation  process  raises  many  issues  which  must  be  understood  in  order  to 
optimize  the  perfOTmance  of  the  model.  The  final  design,  tuning  and  intelligent  use  of  the 
assimilation-fmtcast  system  requires  knowledge  of: 

1)  the  relative  merits  of  different  assimilation  strategies. 

2)  the  cost/benefit  ratio  of  an  extensive  data  assimilation  cycle  versus  a  "cold-stan" 
initialization. 

3)  a  method  to  estimate  the  optimal  value  for  each  of  the  nudging  coefficients. 

4)  the  error  growth  in  MASS  simulations  during  the  forecast  period. 

5)  the  relative  impact  of  each  data  type  on  the  MASS  mesoscale  forecasts. 

An  observing  system  simulation  experiment  (OSSE)  is  an  excellent  tool  for  addressing 
these  issues.  An  OSSE  makes  use  of  a  high  resolution  numerical  model  simulation  to  provide  a 
dynamically  consistent  "surrogate  atmosphere  simulation"  (SAS).  Simulated  observational  data  is 
then  extracted  ftom  the  SAS  and  objectively  analyzed  to  create  a  model  initialization  dataset  and/or 
gr^ed  data  to  be  assimilated  throu^  nudging.  Experimental  simulations  are  then  run  on  a  coarser 
grid  which  is  entirely  contained  within  the  SAS  grid.  The  main  advantage  of  OSSEs  over 
conventional  experimental  simulations  is  that  the  state  of  the  "surrogate  atmosphere”  is  known 
exactly  while  the  state  of  the  real  atmosphere  is  not.  Therefore,  OSSE  simulations  can  be 
rigorously  compared  to  the  "true"  state  of  the  atmosphere.  This  allows  for  the  careful  evaluation  of 
the  performance  of  various  data  assimilation  schemes.  Unfortunately,  there  are  also  several 
disadvantages  to  OSSEs.  For  example,  many  errors  which  find  their  way  into  actual  observations 
cannot  be  included  in  OSSEs.  These  include  errors  introduced  through  the  aliasing  of  small  scale 
features  that  are  not  resolved  on  the  OSSE  grid,  occasional  large  errors  related  to  equipment 
malfunction,  and  systematic  errors  at  one  or  more  observing  locations.  In  addition,  the  nature  of 
data  extraction  in  an  OSSE  makes  it  very  difficult  to  produce  certain  types  of  simulated  data,  such 
as  satellite  and  radar  data  and  visual  cloud  observations.  Finally,  the  "surrogate  atmo^here"  only 
resembles  the  real  atmosphere  to  the  extent  that  the  numerical  model  effectively  simulates  all 
atmospheric  processes.  The  net  result  is  that  OSSE  experimental  simulations  are  expected  to 
parallel  the  SAS  more  closely  than  conventional  simulations  parallel  the  real  atmosphere.  Still,  the 
possibility  of  detailed  comparisons  between  the  experimental  simulations  and  a  well  known 
"verification"  make  an  OSSE  a  worthwhile  endeavor. 

Stauffer  et  al.,  (1991)  point  out  that  precipitation  is  the  best  variable  for  verifying  the 
effects  of  nudging  because  it  is  not  assimilat^  into  the  model  itself  and  it  is  sensitive  to  many 
dynamic  and  thermodynamic  processes  on  all  scales  of  motion.  Therefore,  precipitation  will  be  the 
primary  verification  variable.  However,  since  the  simulations  presented  in  Stauffer  et  al.,  (1991) 
consist  of  only  a  pre-forecast  period,  they  have  no  reason  to  use  any  variable  which  is  directly 
assimilated  for  verification.  Since  the  simulations  presented  here  include  a  short  pre-forecast  period 
followed  by  a  forecast,  other  variables  including  winds,  temperature  and  sea  level  pressure  will 
prove  useful  for  verification  during  the  forecast  period. 

7.7  Experimental  Design 

The  case  of  the  Midwestern  squall  line  of  5-6  September  1992  was  selected  for  the  OSSE. 
A  brief  review  of  the  case  is  presented  in  Section  7.2.  The  period  of  interest  is  the  12  hour  period 
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beginning  at  0000  UTC  6  September.  The  grid  selected  for  the  "surrog^  atmosphere"  emulation 
(SAS)  measures  77  x  70  with  28.6S  km  resolution  and  is  shown  in  Figure  7-la.  The  simulation 
was  initialized  at  0000  UTC  5  September  and  integnued  for  36  hours. 

The  experimental  simulations  were  inidali^  from  the  SAS  at  0000  UTC  6  September  or  at 
1200  UTC  S  September  (see  Table  7-1  and  Figure  7-2).  All  experimental  simulations  were 
integrated  until  1200  UTC  6  September  on  a  41  x  38  grid  with  40  km  resolution.  The  grid  is 
entirely  contakied  within  the  SAS  grid  and  is  shown  in  Figure  7-lb. 

Simulated  data  was  extracted  from  the  SAS  at  standard  observing  station  locations  and 
times.  Locations  of  the  simulated  surface,  rawinsonde  and  profiler  observations  are  shown  in 
Figure  7-3a.  b  and  c.  respectively.  Surface  profiler  statimis  are  assumed  to  report  on  an  hourly 
basis  and  rawinsonde  stations  at  0000  and  1200  UTC.  It  is  assumed  that  any  given  observing 
station  has  a  2.5%  chance  of  not  reporting  at  a  given  time.  In  tnder  to  more  accurately  simulate  real 
observing  systems,  a  randknn,  nora^y  &tributed  error  with  a  standard  deviation  similar  to  that  of 
the  actud  observing  system  was  ad^  to  each  observation.  The  standard  deviation  for  each 
variable  aixl  (^rserving  system  are  riiown  in  Table  7-2. 

Tabu  7-1.  Experimental  Simulations  performed  in  the  OSSE. 


Simulation 

G  (s-l) 

Data  Assimilated 

Cl 

00/6 

Static  initialization 

C2 

00/6 

— 

Same  as  Cl,  but  constant  boundary  conditions 

C3 

12/5 

— 

Same  as  Cl,  but  24  hour  simulation 

SI 

00/6 

.0003 

Surface  u,  v,  T,  q  (XX)0  -  0200  UTC 

S2 

12/5 

.0003 

Surface  u,  v,  T,  q  1900  -  0200  UTC 

U1 

12/5 

.0003 

Rawinsonde  u,  v,  T,  q  2000  -  00(X)  UTC 

U2 

12/5 

.0009 

Rawinsonde  u,  v,  T,  q  2000  -  0000  UTC 

U3 

12/5 

.0001 

Rawinsonde  u,  v,  T,  q  2000  -  0000  UTC 

U4 

12/5  . 

.0003 

Rawinsonde  T,  q  2000  -  0000  UTC 

U5 

12/5 

.0003 

Rawinsonde  u,  v  2000  -  OOOO  UTC 

PI 

12/5 

.0003 

Profiler  u,  v,  derived  T  19(X)  -  0200  UTC 

P2 

12/5 

.0003 

Profiler  u,  v  1900  -  0200  UTC 

P3 

12/5 

.0003 

Profiler  derived  T 1900  -  0200  UTC 

P4 

12/5 

.0003 

same  as  PI,  but  fill  domain  with  profilers 

P5 

12/5 

.0003 

same  as  P4,  but  time-smooth  data 

P6 

12/5 

.0003 

emv-free  wind  observations  at  every  grid  point 

P7 

12/5 

.0003 

same  as  P6  but  nudge  surface  u,  v,  T  and  q  1900-0200  UTC 

Table  7-2.  RMS  errors  for  simulated  observations. 


Surface  winds 

5% 

Rawinsonde  height  (5(X)  mb) 

10  m 

Surface  temperature 

0.3^ 

Rawinsonde  height  (250  mb) 

20  m 

Surface  RH 

0.02 

Rawinsonde  temperature 

0.4®C 

Surface  pressure 

100  Pa. 

Rawinsonde  winds 

12% 

Profiler  winds 

5% 

Rawinsonde  RH 

0.025 

Profiler  heights 

0.5% 
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Figure  7-1  A  depiction  of  the  geographic  area  covered  by  the  (a)  coarse  mesh  and  (b)  fine  mesh 
components  of  the  OSSE. 
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Figure  7~2  Schematic  depicting  OSSE  design.  A  timeline  chart  which  shows  the  types  of  data 
and  a  schedule  for  assimilating  them  into  the  MASS  model  during  the  Observation  System 
Simulation  Experiment  (OSSE). 
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Surfoce  Station  Locations  Rawinsonde  Station  Locations 


Figure  7-3  Locations  of  simulated  (a)  surface,  (b)  rawinsonde  and  (c)  profiler  observations 
used  in  the  OSSE.  Box  encloses  the  41  x  38  inner  grid. 
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The  experimental  simulations  were  initialized  with  the  simulated  rawinsonde  and  surface 
data  only.  The  first  mss  field  for  the  dau  c'^alysis  was  created  from  the  observations  themselves 
rather  than  frx>m  NGM  or  MASS  model  <ridded  data.  This  prevented  the  introduction  of 
infimnaticHi  other  than  that  contained  in  the  simulated  observations.  All  experimental  simulations 
used  the  1200  UTC  S  Septemba  NGM  simulation  for  boundary  conditions. 

Simulated  rawinsonde,  profiler  and  surface  data  were  also  objectively  analyzed  for  the 
purpose  of  nudging.  Gridded  rawinsonde  data  was  created  at  1200  UTC  S  September  and  0000 
UTC  6  September.  Simulated  profiler  winds  were  analyzed  each  hour  between  1900  UTC  5 
September  and  0200  UTC  6  September.  3-D  temperuure  fields  derived  from  profiler  winds 
(Cram,  et  al.  1991)  were  also  available  at  those  times.  Finally,  gridded  surface  data  were  analyzed 
at  the  same  times  as  profiler  data. 

Seventeen  experimental  simulations  were  run  ^d  are  detailed  in  Table  7- 1  and  Hgure  7-2. 
Column  2  in  the  table  refers  to  the  initialization  time  (0000  UTC  6  September  or  1200  UTC  5 
September)  and  G  in  column  3  refers  to  the  assimilation  coefficient  Simulation  Cl  is  the  control 
simulation  with  only  a  static  initialization  at  0(XX)  UTC  6  September.  Simulation  C2  tests  the 
impact  of  the  boundary  conditions  on  the  simulation.  Simulation  C3  is  identical  to  simuladon  Cl 
except  that  it  is  initialized  at  1200  UTC  5  September.  Simulations  SI  and  S2  are  12  and  24  hour 
simulations,  respectively,  which  assimilate  surface  data.  Simulations  Ul,  U4  and  US  assimilate 
various  combinations  of  rawinsonde  data  (u,  v,  T.  q).  Simulations  U2  and  U3  test  the  effect  of 
varying  the  data  assimilation  coefficient.  G  » .0003  implies  an  e-folding  time  of  just  under  1  hour 
for  the  relaxation  of  the  model  state  to  the  "observed"  state.  Simulations  PI,  P2  and  P3  test  the 
assimilation  of  various  combinations  of  profiler  winds  and  derivnl  heights.  Simulations  P4,  PS 
and  P6  are  attempts  to  improve  the  assimilation  of  profiler  data  including,  adding  profiler 
observations  in  tenons  which  are  not  covered  by  the  current  network  (P4),  smtxxhing  the  profiler 
winds  in  time  (PS)  and  creating  a  database  of  error  free  or  "perfect"  profiler  observations  at  each 
model  grid  point  ^6).  Finally,  simulation  P7  assimilates  all  surface  data,  in  addition  to  "perfect" 
profiler  winds  and  temperature.  The  output  from  each  simulation  will  be  verified  against  the  "true" 
values  of  the  variables  from  the  SAS. 

The  performance  of  the  various  experimental  simulations  were  evaluated  for  a  12  hour 
peri^  beginning  at  0000  UTC  6  September.  Because  they  are  not  directly  assimilated, 
precipitation  and  surface  pressure  were  used  to  assess  the  performance  of  each  simulation.  For 
purposes  of  simplicity,  OOOO  UTC  6  September  will  be  referred  to  as  0  hours,  0600  UTC  will  ^ 
referred  to  as  6  hours,  etc. 


7.2  Review  of  the  6  September  1992  case 

The  case  of  the  Midwestern  squall  line  of  S-6  September  1992  was  selected  for  the  OSSE. 
A  squall  line  developed  ahead  of  a  cold  frontal  boundary  during  the  afternoon  of  3  September  and 
propagated  eastward  through  the  eastern  plains  states  before  dissipating  by  about  1200  UTC  6 
September  (Figure  7-4).  At  0035  UTC  (Figure  7-4a)  the  squall  line  extended  from  nonhwest 
Wisconsin  to  west  Texas.  It  was  unbroken  along  its  northern  and  central  sections  while  the 
southern  part  was  broken.  At  0035  and  0635  (Figure  7-4a,  b),  an  unbroken  line  of  precipitation 
extended  from  the  Upper  Peninsula  of  Michigan  to  Oklahoma.  By  0935  (Figure  7-4c)  the  northern 
part  of  the  squall  line  propagated  to  the  western  shore  of  Lake  Michigan  and  began  to  dissipate, 
while  intense  convection  developed  in  eastern  Oklahoma  at  the  southern  end  of  die  line.  At  1235 
(Figure  7-4d)  only  remnants  of  the  squall  line  remained. 
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Figure  7-4  Manually-digitized  radar  plots  for  the  OSSE.  The  times  are  (a)  0335,  (b)  0635,  (c) 
0^5,  and  (d)  1235  UTC  6  September.  Each  level  of  shading  represents  one  VIP  level. 
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73  The  Surrogate  Atmosphere  Simulation 

Figure  7-S  shows  3'hourly  p^pitadon  totals  from  the  SAS  during  the  period  of  interest. 
During  the  12  hour  period  a  squidl  line  propagated  eastward  through  the  ^main  ahead  of  a  cold 
ftont  Initially,  the  ^uall  line  was  locat^  along  an  axis  which  ran  q>proximately  from  extreme 
southwestern  Oklahoma  to  northwestern  Wisconsin.  Over  the  course  of  the  simulation,  the  squall 
line  propagated  eastward  several  hundred  kilometers  while  the  northern  part  of  the  line  dissipated 
after  about  0600  UTC.  A  compari^  of  Figure  7-S  with  Figure  7-4  shows  a  remarkable  agreement 
between  the  SAS  and  Manually  Digitized  Radar  plots  for  the  same  period. 


7.4  Results 

Several  statistical  scores  are  used  to  compare  the  various  simulations.  They  include  the 
0.25  cm  threshold  precipitation  threat  scores  for  each  3  hour  period  and  for  the  entire  12  hour 
period,  and  die  root  mean  square  (RMS)  surface  pressure,  temperature  and  wind  (average  of  u  and 
v)  errcvs  at  0, 6  and  12  hours.  Threat  scores  and  RMS  errors  are  calculated  using  a  subset  of  the 
inner  grid  that  does  not  include  the  four  rows  (or  columns)  of  grid  points  nearest  to  each  boundary. 
The  threat  score  is  determined  by 


TS  = - CEA - 

(FA  +  OA  -  CFA) 


(7-1) 


where  FA  is  the  area  where  the  precipitation  is  forecast  to  equal  or  exceed  the  threshold,  OA  is  the 
area  where  the  observed  precipitation  equals  or  exceeds  the  threshold,  and  CFA  is  the  correctly 
forecast  area  or  the  region  where  both  forecast  and  observed  values  exceed  tire  threshold.  Since  tire 
threat  score  penalizes  a  simulation  for  erroneously  predicting  precipitation  either  above  or  below 
the  threshold  value,  a  perfect  score  of  1.0,  or  even  a  near  perfect  score,  is  nearly  impossible  to 
attain.  The  relatively  shon  3  hour  verification  period  in  this  study,  combined  with  relatively  high 
horizontal  resolution,  increases  the  challenge  funher.  The  3  hour  threat  score  presented  here 
effectively  test  the  nxxlers  skill  at  predicting  meso-B  scale  piecipitaticm  events  in  iq>ace  and  time. 


7.4. 1  Influence  of  boundary  conditions 

Figure  7-6  shows  the  surface  pressure  difference  between  simulations  Cl  and  C2.  At  3 
hours,  differences  are  relatively  small  except  near  the  boundaries.  By  6  hours,  significant 
differences  in  surface  pressure  and  in  surface  pressure  gradient,  as  shown  by  gra^ents  in  the 
difference  field,  are  evident  throughout  most  of  the  model  domain.  These  differences  continue  to 
increase  through  12  hours.  An  examination  of  Table  7-3  shows  that  threat  scores  for  the  two 
simulations  are  similar  through  6  hours.  After  that,  simulation  C2  deteriorates  at  9  hours,  although 
it  outperforms  simulation  Cl  during  the  Hnal  3  hour  period.  Since  the  constant  boundary 
conditions  in  C2  are  certainly  inferior  to  the  boundary  conditions  in  simulation  Cl,  this  can  only  be 
attributed  to  pure  chance.  At  6  hours,  RMS  wind  and  pressure  errors  for  simulation  C2  are 
somewhat  larger  than  those  for  simulation  Cl.  By  12  hours  they  are  significantly  larger. 
Surprisingly,  temperature  errors  are  smaller  than  for  simulation  Cl.  If  any  conclusion  can  be 
drawn  from  this  information,  it  is  that  up  until  6  hours  differences  between  the  experimental 
simulations  are  primarily  due  to  differences  in  data  assimilation.  The  influence  of  the  boundary 
conditions  becomes  significant  after  6  hours  and  continues  to  increase  until  12  hours. 
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Surrogate  Atmosphere  Simulatioii 


ending  09  UTC  ending  12  UTC 


Figure  7~S  Three  hour  accumulated  precipitation  in  hundredths  of  an  inch  ending  at  (a)  0300 
UTC.  (b)  0600  UTC.  (c)  0900  UTC  and  (d)  1200  UTC  06  September  1992 for  the  surrogate 
atmosphere  simulation.  Contour  interval  is  0.10  in.  starting  at  0.01  in. 
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Figure  7-6  Surface  pressure  difference  between  simulations  Cl  and  C2  as  (a)  0300  UTC,  (b) 
0600  UTC,  (c)  0900  UTC  and  (d)  1200  UTC  6  September  1992.  Positive  values  indicate  thcu 
pressure  is  higher  in  simulation  Cl. 
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7.4.2  Statistical  evaluation 

An  examination  of  Table  7-3  shows  that  simulation  SI  is  vinually  identical  to  simulation 
Cl.  The  same  can  be  said  for  simulation  S2  and  simulation  C3.  The  assimilation  of  surface  data 
has  virtually  no  impact  on  the  model's  performance,  at  least  for  the  case  examined  here.  If 
simulation  SI  is  ignored,  simulation  Cl,  the  control  simulation,  has  the  best  overall  threat  scores. 
After  6  hours,  simulation  Cl  outperforms  all  other  experimental  simulations.  Before  6  hours 
however,  several  experiments  which  assimilate  rawinsonde  data  are  superior  to  simulation  Cl. 
These  include  simulation  U1  and  U2.  Simulation  U2  has  the  best  overall  performance  during  this 
period.  Simulations  U3,  U4  and  U5  performed  poorly  relative  to  simulation  Cl  except  during  the 
hrst  period  These  three  simulations  did  produce  a  better  forecast  than  simulation  FC  which  was 
initialized  at  the  same  time  but  did  not  include  any  nudging.  This  indicates  Liat  the  assimilation  of 
data  with  a  weaker  nudging  coefficient  (0.0001)  or  the  assimilation  of  only  wind  or  only  mass 
data,  do  improve  the  simulation  of  precipitation  slightly,  although  not  enough  to  outperform 
simulation  Cl  which  is  initialized  12  hours  later.  Simulations  PI  through  P6  showed  little,  if  any. 
improvement  over  simulation  C3.  Simulation  P7  was  substantially  better  than  simulation  C3, 
aldough  the  precipitadon  was  not  simulated  as  well  as  simuladon  Cl,  U1  or  U2. 


7.4.3  Discussion 

An  examination  of  Table  7-3  and  Figure  7-7  shows  that  simulation  Cl  (the  control 
simuladon)  suffers  from  a  precipitation  spin-up  problem.  Bias  scores  for  the  first  3  hour  period 
(not  shown)  confirm  that  the  total  area  of  precipitadon  is  significandy  smaller  (0.66  for  0. 1  cm 
threshold  and  0.S3  for  0.25  cm)  than  for  the  SAS,  while  for  many  of  the  nudged  simulations  the 
bias  scores  are  significantly  closer  to  1.0.  After  the  first  period,  simuladon  Cl  generally  performs 
as  well  as  or  better  than  any  of  the  other  simulations,  especially  after  6  hours  when  the  threat 
scores  for  simulations  with  upper  air  nudging  drop  to  near  0.  A  subjective  comparison  of  the 
precipitadon  field  from  simulation  Cl  (Figure  7-7)  with  that  from  the  SAS  (Figure  7-5)  shows  that 


Table  7-3.  Threat  scores  and  RMS  errors  for  experimental  simulations. 
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Simulati<ni  Cl 


ending  09  UTC  ending  12  UTC 


Figure  7-7  Three  hour  accumulated  precipitation  in  hundredths  of  an  inch  ending  at  (a)  0300 
UTC,  (b)  0600  UTC,  (c)  0900  UTC  and  (d)  1200  UTC  06  September  1992  for  simulation  Cl. 
Contours  as  in  Figure  7-5. 
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during  the  first  6  hours,  simulation  Cl  fails  to  generate  precipitation  along  the  northern  pan  of  the 
squall  line.  After  6  hours,  the  general  representation  of  the  precipitation  field  is  good,  although 
specific  details  such  as  the  heavy  precipitation  in  southeast  Kansas  between  6  and  9  hours,  or  in 
southeast  Missouri  between  9  and  12  hours,  are  missed  by  simulation  Cl. 

Simulation  SI,  which  assimilates  surface  data,  is  nearly  identical  to  simulation  Cl.  The 
same  is  true  for  simula^ons  S2  uul  C3.  This  implies  that  in  this  case  and  others  where  the  majority 
of  the  forcing  is  on  the  synoptic  scale,  siuface  nudging  alone  has  litde  or  no  beneficial  effect 

Unlike  surface  nudging,  upper  air  nudging  (simulations  U1  through  US)  has  a  considerable 
effect  on  die  simulations.  The  precipitation  field  from  simulation  U2  is  shown  in  Figure  7-8.  Up 
until  6  hours,  simulation  U2  does  a  remarkable  job  of  simulating  the  precipitation  from  the  SAS. 
Simulation  U1  (not  shown)  does  nearly  as  well.  The  stronger  nudging  coefficient  in  simulation  U2 
results  in  a  less  bdanced  state  at  0  hours  as  evidenced  by  larger  RMS  picture  and  wind  enors.  By 
6  hours  though,  errors  for  the  two  simulations  are  similar.  Still,  it  is  not  inconceivable  that  under 
some  circumstances,  the  stronger  nudging  coefficient  in  simulation  U2  may  cause  problems  by 
pushing  the  simulation  too  far  out  of  balance.  Threat  scores  for  simulation  U3  in^cate  that  a 
nudging  coefficient  of  0.0001  is  too  small  to  adequately  introduce  rawinsonde  data  into  the 
simulation.  Simulation  U4  which  assimilates  temperature  and  moisture  only  and  simulation  US, 
which  assimilates  winds  only,  both  perform  rather  poorly.  It  seems  that,  at  least  for  the  present 
case,  both  the  wind  and  mass  fields  need  to  be  assimilated  to  improve  the  simulation. 

After  6  hours,  simulations  U1  and  U2  weaken  the  squall  line  too  quickly  and  develop 
phase  errors.  In  fact,  simulation  Cl  outperforms  all  nudged  simulations  after  6  hours.  In  fact, 
RMS  errors  for  simulations  U1  through  U5  are  larger  throughout  the  period  of  interest,  except  in 
the  case  of  surface  pressure  at  6  hours.  It  seems  that  rawinsonde  nudging  improves  the  simulation 
during  the  first  6  hours  by  eliminating  the  spin-up  problem.  Simulation  U1  also  provides  a  more 
balanced  initial  state  as  shown  by  the  fact  that  the  RMS  surface  pressure  error  increases  by  only 
0.13  mb  in  the  first  3  hours  compared  to  0.50  mb  for  simulation  Cl.  After  6  hours,  larger  initi^ 
RMS  wind  and  mass  field  errors  contribute  to  a  degradation  of  simulations  U1  and  U2. 

The  assimilation  of  simulated  profiler  data  produced  disappointing  results.  The  squall  line 
lagged  several  hundred  kilometers  behind  its  "observed”  position  in  simulation  PI  (not  shown). 
Tli^t  scores  are  poor  throughout  the  simulation.  RMS  errors  are  larger  than  for  simulation  U1 
except  for  winds  at  0  hours. 

Two  possible  explanations  for  the  poor  performance  are;  (1)  the  profiler  network  extends 
over  only  a  portion  of  the  model  domain  (Figure  7-2c).  The  only  result  of  profiler  nudging  for  the 
regions  outside  the  profiler  network  is  that  there  are  12  addition^  hours  during  which  error  growth 
can  occur.  (2)  the  profiler  network  suffers  from  an  inconsistency  between  low  spatial  and  high 
temporal  resolution.  This  problem  is  exacerbated  because  the  simulated  profiler  data  was  created 
from  a  single  observation  time,  unlike  the  actual  profiler  network  which  averages  5  minute  data  to 
produce  hourly  observations.  The  result  is  a  temporally  noisy  analysis  which  is  then  nudged  into 
the  simulation. 

Several  experimental  simulations  were  run  to  test  these  possibilities.  In  simulation  P4, 
additional  simulate  profiler  stations  were  added  so  that  the  entire  model  domain  was  adequately 
covered  by  the  data.  Simulation  PS  was  identical  to  simulation  P4,  but  a  time  smoother  was  applied 
to  the  dau  to  remove  features  which  occurred  on  a  time  scale  of  less  than  three  hours.  Neither 
produced  any  significant  improvement  over  simulation  PI  (see  Table  7-3).  As  a  final  test, 
simulated  profiler  data  was  extracted  at  each  model  gridpoint  (simulation  P6).  For  this  simulation, 
the  RMS  errors  were  set  to  0  (see  Table  7-2),  and  the  time  tendency  term  was  not  ignored  in  the 
derivation  of  temperatures  from  the  wind  field.  Even  with  "perfect "  profiler  data  at  every  model 
gridpoint,  there  was  no  significant  improvement  on  simulation  PI.  Simulation  P7  was  identical  to 
simulation  P6  with  the  addition  of  nudging  to  surface  wind,  temperature  and  moisture.  This 
simulation  was  the  only  one  for  which  profiler  nudging  produced  a  squall  line  in  approximately  the 
same  location  as  for  the  SAS.  Even  so,  the  statistics  for  simulation  P7  are  considerably  worse  than 
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Simulatiini  U2 


ending  09  UTC  ending  12  UTC 


Figure  7-8  Three  hour  accumulated  precipitation  in  hundredths  of  an  inch  ending  at  (a)  0300 
UTC,  (b)  0600  UTC.  (c)  0900  UTC  and  (d)  1200  UTC  06  September  1992 for  the  simulation  U2. 
Contours  as  in  Figure  7-5. 
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those  for  simulation  Ul. 

Given  that  the  simulated  profiler  wind  field  was  essentially  perfect,  the  one  remaining 
explanation  for  the  poor  performance  of  profiler  data  assimilation  relative  to  rawinsonde  data 
assimilation  is  the  method  of  deriving  the  temperature  field  from  the  wind  field.  Figure  7>9  shows 
the  700  mb  temperature  analyses  from  the  simulated  rawinsonde  temperature  observations,  the 
temperature  fields  derived  from  the  simulated  “standard"  profiler  wind  observations  us^  in 
simulation  P4  and  fiom  the  "perfect"  profiler  wind  observations  used  in  simulations  P6  and  P7. 
The  figure  also  shows  the  "verification"  from  the  SAS.  At  700  mb,  all  of  the  analyses 
underestimate  the  temperature  in  the  warm  tongue  which  extends  from  Oklahoma  into  Wisconsin. 
The  standard  profiler  analysis  underestimates  the  temperature  the  most  followed  by  the  perfect 
profiler  analysis  and  the  rawinsonde  analysis.  The  temperature  field  derived  from  the  "perfect" 
profiler  anal^^s  resolves  the  small  scale  feanire  from  the  SAS  very  well,  although  the  tenqieranires 
are  ever^here  about  2  K  colder  than  the  SAS.  The  temperature  field  from  the  standard  profiler 
observations  resolves  the  details  of  the  SAS  temperanire  field  substantially  better  than  the 
rawinsonde  analysis,  but  does  not  show  some  of  the  smaller  scale  features  evident  in  the  perfect 
analysis.  It  too  is  generally  several  degrees  too  cold.  The  rawinsonde  analysis  is  the  smoodiest  of 
the  three;  however,  it  underestimates  temperatures  the  least,  especially  in  the  tongue  of  warm  air 
which  extends  from  Oklahoma  to  Wisconsin.  At  350  mb  (not  shown),  the  profiler  ansdyses  tend  to 
overestimate  the  temperature  relative  to  the  SAS  while  the  rawinsonde  analysis,  although  it  doesn't 
reveal  as  many  details,  shows  no  bias  in  the  temperature  field.  The  net  result  of  the  bias  introduced 
into  the  temperature  field  by  nudging  to  derived  profiler  temperature  field  is  that  the  lapse  rate  is 
reduced  in  and  around  the  squall  Une  (Figure  7-10)  which  may  explain  the  poor  performance  when 
profiler  data  is  assimilated. 

What  is  the  reason  for  the  excellent  representation  of  small  scale  features  in  the  derived 
profiler  temperature  field,  but  the  large  overall  bias?  When  the  temperature  field  is  derived  firom  the 
profiler  wind  fields,  the  technique  of  Cram  et  al.  (1991)  actually  calculates  the  height  field.  Two 
possible  sources  of  error  are  0)  the  source  of  the  height  field  used  for  boundary  conditions  and  (2) 
the  amplification  of  small  height  errors  when  the  temperature  field  is  calculated  from  the  thickness 
equation.  In  the  derivation  of  temperatures  from  the  perfect  profiler  data,  model  output  from  the 
SAS  was  used  for  bound^  conditions.  This  rules  out  (1).  (2),  however,  deserves  some 
consideration.  Three  dimensional  nudging  data  was  analyzed  on  pressure  surfaces  spaced  every  50 
mb  (about  500-1000  m)  in  the  vertiesd.  Therefore,  a  5  to  10  meter  error  in  the  height  field  at  any 
one  pressure  level  would  produce  about  a  1%  error  in  the  thickness  between  analysis  levels  or 
about  a  2  to  3  K  error  in  the  temperature  field!  It  appears  that  the  technique  for  deriving  the 
temperature  field  from  the  wind  field  produces  a  temperature  analysis  that  is  considerably  poorer 
than  the  rawinsonde  temperature  analysis.  Perhaps  there  is  a  way  to  merge  profiler  wind  data  with 
rawinsonde  temperature  data,  although  initial  attempts  at  this  have  not  pr^uced  good  results. 

The  addition  of  surface  nudging  is  all  that  distinguishes  simulation  P7  from  simulation  P6, 
but  there  is  a  significant  difference  between  the  two  simulations.  Sirhulation  P7  is  the  only  "P" 
series  simulation  in  which  the  squall  line  is  well  represented  during  any  part  of  the  verification 
period  (not  shown).  One  possible  explanation  is  that  the  nudging  of  surface  temperature 
observations  into  the  PBL  helps  to  overcome  the  tendency  of  profiler  temperature  nudging  to 
stabilize  the  atmosphere.  When  surface  temperature  nudging  is  added,  all  gridpoints  within  the 
model-defined  PBL  are  nudged  to  surface  data  rather  than  to  the  "perfect"  profiler  data.  Of  course, 
wind  or  moisture  nudging  in  the  PBL  may  have  also  played  a  role.  In  any  case,  the  combination  of 
surface  nudging  and  rawinsonde  nudging  (not  shown)  produces  nearly  identical  results  to  that  of 
rawinsonde  nudging  alone.  So  it  appears  that  the  surface  data  is  helpful  only  because  it 
compensates  for  a  weakness  in  the  profiler  nudging  and  not  because  it  contains  any  information 
which  enhances  a  situation  in  which  the  upper  level  dynamics  are  well  represented. 
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Figure  7-9  700  mb  temperature  (K)  valid  at  0000  UTC  06  September  1992  (a)  from  SAS,  (b) 
from  simulated  rawinsonde  analysis,  (c)  derived  from  simulated  wind  profiler  analysis  and  (d) 
derived  from  sinudated  "perfect"  wind  profiler  analysis.  Contour  interval  is  1 K. 
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Figure  7-10  Vertical  cross  section  of  potential  temperature  (K)  at  0000  UTC  6  September  1992 
along  the  path  indicated  in  (a).  Cross-sections  are  from  (b)  the  SAS,  (c)  simulation  P6.  and  (d) 
simulation  Ul.  The  left  edge  of  (b),  (c)  and  (d)  correspond  to  the  southeast  edge  of  the  path  in  (a). 
Contour  interval  is  3  K. 
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7.5  Conclusions 

The  results  of  the  OSSE  indicate  that  the  main  benefit  of  nudging  is  in  the  reduction  of  the 
spin-up  pn^em  during  the  initial  3  to  6  hours  of  a  simulation.  After  this  time,  nudged  siinulations 
are  infenor  to  a  stadcally  initialized  simulation  because  the  nudged  simulation  must  be  initialized  12 
hours  earlier,  at  least  with  the  current  availability- of  complete  3-diinensional  observational  datasets. 
The  results  of  the  case  of  S-6  September  1992  indicate  that  the  assimilation  of  rawinsonde  data 
pr^uces  the  greatest  improvement,  while  surface  data  is  less  valuable  and  profiler  data  suffers 
from  a  lack  of  information  about  the  mass  field.  Since  the  OSSE  consisted  of  only  1  case,  these 
results  are  not  necessarily  true  for  all  conceivable  situations.  For  cases  with  we^  forcing,  the 
assimilation  of  rawinsonde  d***  may  be  less  impcntant  while  surface  data  assimilation  may  iuve  a 
greater  impact  Stauffer  et  al.  (1991)  found  this  to  be  true  in  their  weak  synoptic  fmcing  case. 

An  impOTtant  consideration  in  the  assimilation  of  surface  data  is  that  the  state  of  the  PEL, 
especially  temperature,  respcmds  to  external  forcing  on  very  shon  time  scales.  For  this  reason,  if 
the  farcing  fiinctkms  in  the  PEL  are  not  properly  simulated,  the  assimilation  of  surface  data  into  the 
PEL  may  improve  the  surface  representadveness  of  the  model  PEL  only  undl  the  nudging  period 
ends.  It  may  be  more  valuable  to  concentrate  on  accurately  representing  atmospheric 
transmissivity,  soil  tix>isture  and  other  variables  which  determine  the  surface  forcing  function.  The 
assimilation  of  Manually  Digitized  Radar  (MDR)  data  and  the  synthetic  relative  humidity  scheme 
(see  Sections  6.4  and  6.5),  which  unfortunately  could  not  be  tested  in  an  OSSE,  may  help  to 
improve  the  model's  representation  of  cloud  cover  which  implies  a  better  representation  of  the 
radiative  component  of  the  surface  forcing  function.  It  may  dso  be  useful  to  adjust  the  vertical 
structure  of  the  assimilation  of  surface  data  into  the  PEL  to  the  particular  PEL  structure  (free 
convection,  stable,  etc.)  assigned  by  the  Elackadar  PEL  scheme  at  each  gridpoint.  This  would  be 
consistent  with  Stauffer  et  al.,  (1991)  who  emphasized  that  the  method  of  assimilation  of  surface 
data  must  not  conflict  witii  the  model  PEL  physics. 

The  assimilation  of  profiler  data  could  prove  to  be  quite  valuable  if  a  method  of  removing 
the  temperature  bias  in  the  derived  temperature  field  could  be  derived.  One  such  method  is  to  use 
rawinsonde  temperature  observations  to  anchor  the  profiler  derived  temperature  field. 
Unfortunately,  rawinsonde  data  is  typically  available  only  every  12  hours  compart  to  every  hour 
for  profiler  data.  Still,  if  the  temperature  biases  are  relatively  constant  with  time,  it  may  be  possible 
to  remove  the  temperature  bias  from  the  1200  and  0000  UTC  profiler  observations  and  then  apply 
the  same  corrections  to  other  times.  If  this  technique  is  not  successful,  it  will  be  necessary  to  wait 
until  3-dimensional  temperature  data  is  available  at  the  same  frequency  as  profiler  wind  da^ 

Given  an  abundance  of  time  an  other  resources,  OSSEs  have  the  potential  to  expand 
understanding  of  data  assimilation  much  more  than  the  limited  experiment  presented  here.  Several 
potentially  valuable  experiments  must  wait  for  another  opportunity.  They  include; 

1)  The  use  of  e  ations  from  the  SAS,  rather  than  output  from  a  numerical 
model,  to  provide  boundary  condition  information.  This  may  help  to  reduce  the 
degradation  of  the  simulation  which  occurs  when  forecast  boundary  condition 
da^  propagates  into  the  domain  during  the  data  assimilation  period. 

2)  Creating  a  technique  to  simulate  MDR.  IR  satellite  and  other  non-standard  data 
frixn  the  SAS  so  that  the  MDR  assimilation  scheme  and  the  synthetic  RH  scheme 
could  be  tested. 

3)  Testing  the  effects  of  temporal  and  spatial  distribution  of  the  simulated 
observations  on  the  simulation. 
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8.  Results  of  36-Case  Evaluation 

Many  significant  changes  to  the  modd  physics  have  been  made  in  the  course  of  this 
project  In  onter  to  evaluate  the  state  of  the  entire  MASS  system  at  the  end  of  the  project  a  large 
sample  of  mesoscale  simulations  were  carried  out.  The  simulations  have  been  evaluated  both 
objectively  (Section  82),  and  on  a  subjective  basis  (Section  8.3). 

8.1  Design  of  Evaluation  Sample 

Since  the  objective  was  to  evaluate  the  preprocessor  and  model  in  general,  a  set  of 
experiments  was  designed  to  examine  the  model  using  a  broad  sample  of  cases.  There  were  two 
factors  which  were  varied  to  achieve  the  desir^  diversity;  time  of  year  and  geography.  First, 
three  sets  of  dates  were  selected  fcv  simulation.  Table  8-1  lists  the  dates  for  the  36-case 
simulaticms.  One  winter  and  two  summer  cases  were  chosen  to  rest  the  model  under  different 
meteorological  regimes.  The  criteria  for  choosing  these  dares  were  that  MESO  possessed 
relatively  complete  sets  of  raw  data,  and  that  mereorologically-inreresting  events  occurred.  For 
each  of  the  nine  dates,  simulations  were  conducted  over  four  different  regions  of  the  country. 
Each  simulation  involved  a  24  hr  large  scale  run  with  a  grid  spacing  of  50  km,  beginning  at  (XXX) 
UTC  on  the  day  of  simulation,  followed  by  a  12  hr  nested  simulation  with  a  grid  spacing  of  IS 
km,  beginning  at  1200  UTC.  Table  8-1  lists  the  names  assigned  to  the  various  model  domains, 
and  Figures  8-1  through  8-4  show  these  grids.  In  order  to  closely  examine  the  performance  of 
the  m<^el  (especially  various  aspects  of  the  surface  parameterization),  five  diagnostic  points 
were  chosen  for  each  region  (also  listed  in  Table  8-1  and  shown  in  Figure  8-S).  At  each  of  these 
points,  the  model's  SRPH  scheme  diagnostics  were  turned  on  so  ^at  a  large  set  of  surface 
information  was  printed  to  files  as  the  runs  progressed.  This  data  is  used  in  the  discussions  of 
individual  cases  Mow. 
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TahU  8‘I  Characteristics  of  the  36-Case  Sample 


Dates  (aU  1992 


Aunst22.23.24 


December  10. 11. 12 


Laree  Scale 


Key  Mete(»oloeicai  Event 


Hurricane  Lester  m  Southwest 


MCC  in  Iowa 


Bis  Northeastern  Snowstorm 


ames 


Nest 


Northeast 

New  York 

Southeast 

Gulf 

Southwest 

New  Mexico 

Great  Plains 

Illinois 

Domain 


Locations 


Northeast-New  York 

Albany.  NY  (ALB) 

Bradford.  PA  (BFD) 

Boston.  MA  C^S) 

New  York.  NY  (LGA) 

Rochester.  NY  (ROC) 

Southeast-Gulf 

Greenwood.  MS  (GWO) 

Macomb,  MS  (MCB) 

Montgomery.  AL  (MGM) 

New  Orieans.  LA  (MSY) 

Pensacob.  FL  (PNS) 

Southwest-New  Mexico 

Carisbad,NM(CNM) 

ElPaso.TX(ELP) 

Holloman  AFB.  NM  (HMN) 

Los  Cruces.  NM  (LRU) 

Truth-or-Conseouences.  NM  (TCS) 

Great  Plains-Illinois 

Columbia.  MO  (COU) 

Chicago,  IL  (MDW) 

Evansville,  IN  (EVV) 

Moline,  IL  (MLI) 

West  L^avette.  IN  (LAF) 
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Figure  8-1  (a)  Large  scale  and  (b)  nested  grid  domcdns  for  the  Northeast  and  New  York  grids. 
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Figure  8-2  (a)  Large  scale  and  (b)  nested  grid  domains  for  the  Southeast  and  Gidf  grids. 


SBIR  Phase  0  Final  Rqion 


DAAD07-90-C-0134 


Page  74 


Figure  8-3  (a)  Large  scale  and  (b)  nested  grid  domains  for  the  Southwest  and  New  Mexico 
grids. 
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Figure  8-4  (a)  Large  scale  and  (b)  nested  grid  domains  for  the  Great  Plains  and  Illinois  grids. 
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Figure  8-5  Location  of  diagnostic  points  for  (a)  Northeast-New  York,  grid,  (b)  Southeast-Gulf 
gfid,  (c)  Southwest-New  Mexico  grid,  and  (d)  Great  Plains-Ulinois  grid. 
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8.2  Statistical  Evaluation 

In  order  to  provide  an  objective  picture  of  the  performance  of  the  entire  MASS 
preprocessor  and  modeling  system,  a  statistical  analysis  of  the  36-case  sample  was  conducted. 
The  parameters  used  for  evaluation  are  the  root-mean-Spquare  (RMS)  and  bias  errors  of  various 
meteorological  variables  (pressure,  temperature,  dew  point,  wind  speol  and  wind  direction),  and 
precipitation  THREAT  scores.  The  THUIEAT  score  is  computed  by  comparing  observed  station 
precipitation  with  model  precipitation  at  nearby  points.  It  is  given  as 

'  6-^(f-CT)  ■ 

where  CF  is  the  number  of  station  locations  (THREAT  scores  are  often  calculated  from  areas  or 
numbers  of  trKxlel  grid  points  rather  than  station  locations)  at  which  the  model  has  correctly 
forecasted  a  given  amount  of  precipitation  over  a  specified  time  period,  O  is  the  observed 
number  of  stations  which  meet  the  same  criteria,  and  F  is  the  total  number  of  stations  which  the 
model  predicts  to  meet  the  criteria  (both  correctly  and  incorrectly).  The  THREAT  score  both 
rewards  correct  forecasts  and  penalizes  incorrect  forecasts  (the  term  in  parentheses  in  the 
denominator),  so  it  is  a  stringent  test  of  forecast  ability.  A  THREAT  score  of  one  would 
represent  a  perfect  forecast. 

The  data  was  compiled  in  various  ways  in  an  attempt  to  give  insight  into  the 
characteristics  of  the  model.  First,  the  statistical  data  from  each  of  the  case  studies  was  averaged 
for  each  of  the  eight  different  grids  (four  large  scale  and  four  nested).  Figure  8-6  shows  the  RMS 
and  BIAS  errors  for  pressure  for  each  of  the  eight  grids.  The  BIAS  plot  shows  that  the  pressure 
is  systematically  decreasing  on  every  grid.  This  represents  a  loss  of  mass  over  the  entire  grid;  it 
is  a  characteristic  problem  in  limited  area  models.  The  problem  is  especially  significant  for  the 
nested  grids,  because  the  mass  flux  across  the  domain  boundaries  is  very  sensitive  to  the 
formulation  of  the  lateral  boundaiy  conditions.  Figure  8-7  shows  the  temperature  characteristics 
of  the  model.  The  Southwest  and  New  Mexico  grids  have  the  largest  errors,  and  also  tend  to  be 
too  cold.  The  probable  reason  for  this  is  that  in  mountainous  tenrain,  meteorological  stations  tend 
to  be  located  in  population  centers  which  are  usually  at  lower  elevations  than  the  average  of  the 
local  terrain  (at  the  foot  of  mountains,  in  river  valleys).  The  model  terrain  is  averaged  over  an 
area  large  enough  that  the  nearest  model  point  to  a  station  location  is  almost  always  higher  than 
the  station,  and  is  normally  cooler  as  a  result.  In  addition,  the  lateral  diffusion  of  temperature  (to 
maintain  numerical  stability)  results  in  some  systematic  temperature  errors  in  complex  terrain. 
Figure  8-8  shows  the  average  dew  point  errors.  The  model  clearly  has  a  moist  bias,  which  is 
probably  due  to  a  systematic  overestimation  of  evapotranspiration  in  the  model's  surface  energy 
budget,  especially  in  the  winter  simulations.  More  work  is  needed  to  make  the 
evapotranspiration  scheme  more  accurate.  The  dry  bias  on  the  New  Mexico  grid  probably  is 
caused  by  a  too-dry  soil  moisture  assumption  on  the  grid  (0.1  compared  to  0.2  for  the  other 
grids).  Figure  8-9  shows  the  wind  speed  characteristics.  The  low  errors  on  the  Great  Plains  ^d 
is  probably  due  to  the  fact  that  the  surface  is  quite  uniform  for  that  grid,  with  very  little  coastline 
and  no  terrain  to  produce  complicated  local  wind  circulations. 
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Figure  8-6  RMS  and  BIAS  errors  for  pressure  (mb)  averaged  over  all  of  the  simulation  dates. 
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Figure  8-7  RMS  and  BIAS  errors  for  temperature  (T)  averaged  over  all  of  the  simulanon  dates. 
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Figure  8-^  RMS  and  BIAS  errors  for  dew  point  (  T)  averaged  over  all  of  the  simulation  dates. 
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Figure  8-9  RMS  and  BIAS  errors  for  wind  speed  (kts)  averaged  over  all  of  the  simulation  dates. 
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Hgure  8-10  shows  the  0.01  inch  THREAT  scores  for  each  of  the  large  scale  grids  (many 
of  the  nested  grids  had  completely  dry  simulations,  which  made  it  difficult  to  compute 
meaning^  precipitation  statisdcs).  The  THREAT  scores  are  good  in  general,  although  there  are 
wide  variations  between  simulations,  and  at  various  times  in  the  sunuladons.  Pr^pitadon 
remains  the  most  difficult  meteorological  variable  to  predict  There  were  boNth  very  good  and 
very  poor  precipitation  forecasts  in  the  36-case  sample.  The  Northeast  snowstorm  was  mostly 
well-predicted,  but  the  mesoscale  convective  complex  (MCX!)  over  Iowa  in  the  September  series 
of  runs  was  pooriy  handles,  even  when  the  nested  grid  was  shifted  to  be  centered  over  the  stcmn. 
The  main  reason  was  that  the  forcing  for  the  MCC  was  not  well-resolved  by  conventional 
observations,  even  when  supplemented  by  surfa^  and  radar  nudging,  while  the  synoptic  scale 
forcing  responsible  fOT  the  Northeast  storm  was  well-captured  by  the  inidai  and  boundary 
conditions. 
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Figure  8-10  Precipitation  THREAT  scores  (.01  inch)  for  each  of  the  large  scale  grids,  averaged 
over  all  the  simulation  dates. 

Second,  the  statistical  data  was  divided  by  time  in  the  simulation,  to  see  the  evolution  of 
errors  in  a  set  of  similar  simulations.  The  August  24  case  was  selected  for  analysis,  as  a 
summenime  example  with  a  relatively  complete  dataset.  Figure  8-11  shows  the  development 
with  time  of  the  temperature  RMS  and  BIAS  errors  for  the  large  scale  grids  on  August  24.  The 
temperature  errors  grow  with  time,  and  the  BIAS  has  an  interesting  evolution.  The  temperatures 
are  systematically  warm  in  the  early  morning  hours,  peaking  at  about  the  expected  time  of  the 
morning  temperature  minimums,  while  the  model  is  too  cool  by  several  degrees  in  the  afternoon. 
At  least  a  part  of  the  reason  for  this  is  the  fact  that  the  model’s  first  layer  is  higher  than  the  level 
where  standard  observations  are  taken  (about  7  or  8  m  vs.  2  m  for  observations).  This  means  that 
level  1  in  the  model  does  not  see  the  more  extreme  values  observed  during  strong  radiative 
heating  or  cooling.  The  excess  of  evaporation  during  the  afternoon  may  also  be  a  part  of  the 
explanation. 
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FiputS^ll  Temperature  RMS  and  BIAS  errors  for  large  scale  simuiations  on  August  24, 1992. 
Figure  8-12  shows  the  dew  point  enors.  and  Figure  8-13  shows  wind  qieed  enxxs.  The 

wuid  speed  biases  show  that  there  is  a  pronounced  tendncy  in  the  model  for  winds  which  are  too 
strong  *t  night  and  stightly  too  weak  during  the  afternoon.  The  above  riispiigjftnn  of  temperature 
OTors  is  also  relevant  to  the  speed.  The  first  layer  needs  to  become  more  decoupled  from 
me  TOt  of  the  atmosphere  at  ni^t,  so  dot  radiative  cooling  can  drive  the  temperature  minimum  a 
few  degrees  lower. 
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Figun  8-12  Dew  point  RMS  and  BIAS  errors  for  large  scale  simulations  on  August 24, 1992. 
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Figure  8-13  Wind  speed  RMS  and  BIAS  errors  for  large  scale  simulations  on  August  24. 1992. 

A  close  examination  of  the  December  1 1  winter  case  revealed  some  interesting  features. 
The  evolution  of  temperature  errors  (Figure  8-14)  is  quite  different  than  in  the  summer  case 
(Figure  8-11).  The  temperature  errors  are  smaller  in  general,  because  the  situation  is  more 
controlled  by  large  scale  processes  such  as  advection,  than  by  local  boundary  layer  processes. 
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hours  after  initialization 


Figure  8‘t4  Tempentture  RMS  and  BIAS  errors  for  the  two  northern  large  scale  simulations 
(Nmtheast  and  Great  PUuns)  on  December  II,  1992. 

Figure  8-15  shows  the  dew  point  errors  for  the  same  December  11  case.  The 
overesdmation  of  evaporation  is  especially  troublesome  in  the  winter,  when  very  little 
evaporation  should  be  tiding  place. 


I 

I 


RMS 
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hours  after  initialization 

Figure  8-15  Dew  point  RMS  and  BIAS  errors  for  the  two  northern  large  scale  simulations 
(Northeast  and  Great  Plains)  on  December  11, 1^2. 


83  Simulation  Examples 


8.3.1  December  1 1  Northeast  Case:  The  First  "Storm  of  the  Century" 

A  major  storm  which  was  called  by  some  the  "Storm  of  the  Century"  hit  the  northeastern 
U.S.  over  the  period  from  December  10  to  December  12,  1992.  “This  storm  was  later 
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overshadowed  by  an  even  stronger  storm  which  battered  the  entire  East  Coast  in  March,  1993; 
this  storm  was  also  dubbed  the  Strain  of  the  Century,  with  somewhat  mrae  justification. 

In  the  December  strain,  a  strong  surface  low  pressure  center  just  off  the  East  Coast 
associated  with  a  vigraous  upper  level  baroclinic  wave  produced  heavy  rainfall  and  coastal 
flooding  in  New  Jersey  and  New  York  City,  and  heavy  snowfall  inland,  especially  in 
Pennsylvania,  New  York  State  and  Massachusetts.  In  addition  to  the  cases  for  this  same  period 
which  were  a  part  of  the  36-case  sample,  a  special  set  of  MASS  simulations  was  praformed  in 
order  to  focus  the  nested  simulation  on  a  different  twelve  hour  period.  A  40  km  simulation 
began  at  1200  UTC  10  December  and  ran  for  36  hr.  A  nested  IS  km  simulation  used  lateral 
boundary  conditions  from  the  40  km  run.  It  began  at  0000  UTC  1 1  December  and  ran  for  IS  hr. 
These  simulations  were  interesting  because  of  both  the  size  and  strength  of  the  storm,  and 
because  of  some  mesoscale  features  of  the  storm. 

Hgure  8-16  shows  the  modeled  evolution  of  the  surface  pressure  field  during  the  40  km 
simulation.  Multiple  areas  of  low  pressure  which  were  initially  large  and  diffuse  merged  into 
one  much  deeper  (991  mb)  low  off  the  coast  of  Virginia  by  about  0900  UTC  1 1  Deceml^.  Hie 
actual  position  of  the  low  was  a  little  further  nonh.  The  evolution  of  model-produced 
precipitation  is  shown  in  Figure  8-17.  The  heavy  precipitation  to  the  north  of  the  low  pressure 
center  was  observed,  although  the  precipitation  pattern  is  also  shifted  somewhat  to  the  south  of 
were  it  was  observed.  Central  and  northwestern  Pennsylvania  received  heavy  snowfall  through 
the  middle  pan  of  the  simulation.  Pittsburgh,  Johnstown,  Erie,  Bradford,  Altoona,  and  State 
College  all  received  very  heavy  snowfall  during  this  period.  The  precipitation  maximum  north 
of  Lake  Ontario  was  also  verified;  the  area  around  Toronto,  ON  experienced  heavy  snowfall. 

Precipitation  distributions  from  the  15  km  simulation  are  shown  in  Figure  8-18.  The 
precipitation  pattern  is  strongly  related  to  the  local  terrain,  with  a  strong  maximum  developing 
over  the  Catskill  Mountains  in  New  York  and  extending  southwestward  into  elevated  areas  in 
northeastern  Pennsylvania.  A  distinct  minimum  occurs  along  the  Hudson  River  valley,  which 
was  observed  to  occur  as  a  result  of  downsloping  flow  down  the  lee  side  of  the  Berkshire 
Mountains  in  western  Massachusetts.  Figure  8-19  demonstrates  that  the  downslopc  flow  as 

well  simulated  by  the  MASS  model.  The  vertical  velocity  field  (to  =  dp/dt)  in  Figure  8- 19b 

clearly  shows  that  the  easterly  winds  have  induced  strong  upward  motion  (negative  to)  on  the 

windward  sides  of  the  Berkshires  and  the  Catskills  and  subsidence  (positive  to)  on  the  lee  sides. 
This  interpretation  is  supported  by  rain  and  snowfall  data  from  hydrological  observations,  which 
are  plotted  in  Figure  8-20.  The  lack  of  snow  in  the  Hudson  Valley  contrasts  with  the  observation 
of  15  inches  of  new  snow  in  the  Catskills.  Table  8-2  lists  12  hr  precipitation  amounts  for 
selected  standard  surface  stations.  The  15  km  simulation  properly  simulated  the  heavy 
precipitation  in  the  New  York  City  area,  and  the  lack  of  precipitation  in  Albany  and  Glens  Falls. 
Neither  simulation  produced  rainfall  as  heavy  as  the  greater  than  two  inch  amounts  reported  in 
New  Jersey  and  eastern  Pennsylvania. 
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TabU  B<2  Observed  and  simulated  precipitation  (inches)  for  selected  cities  fr<m  (XX)0 
me  to  1200  me  II  December  1992 
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Figure  8-16  Mean  sea  level  pressure  (mb)  from  MASS  40  km  simulation  at  (a)  I2(X)  UTC  and 
(b)  1800  UTC  10  December,  (c)  0000  UTC,  (d)  0600  UTC.  and  (e)  1200  UTC  11  December. 
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Figure  8-17  Precipitation  ( inches  of  liquid  water)  from  the  40  km  MASS  simulation  for  3  hr 
periods  ending  at  (a)  1800  UTC  10  December,  (b)  0000  UTC,  (c)  0600  UTC,  and  (d)  1200  UTC 
11  December. 
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Figure  8-18  Precipitation  (inches  of  liquid  water)  from  the  15  km  MASS  simulation  for  3  hr 
periods  ending  at  (a)  0300  UTC,  (b)  0600  UTC,  (c)  0900  UTC,  and  (d)  1200  UTC  11  December. 
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Figure  8-19  (a)  Winds  in  the  lowest  model  layer  (contoured  in  knots)  and  (b)  vertical  velocity 
(to  =  dp/dt)  from  the  MASS  15  km  simulation  at  1200  UTC  11  December  1992.  The  vertical 
velocity  contours  are  in  Pa  s~f'  the  solid  contours  are  positive  w  (downward  motion),  the  dashed 
contours  represent  negative  to  (upward  motion). 
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8.3.2  September  14  Southeast  Case:  Subtle  Mesoscale  Features 

The  day  of  September  14, 1992  was  a  fairly  trai^uil  day  in  the  Southeast,  with  an  upper 
level  ridge  of  high  pressure  dominating  the  synopdc  situation.  The  set  of  MASS  simulations 
which  were  carried  out  on  this  day  contains  some  subtle  features  which  were  quite  well- 
simulated. 

Figure  8-21  shows  model  fields  halfway  through  the  nested  MASS  simulation  at  1800 
UTC  14  September.  Three  mesoscale  features  are  interesting: 

(1)  The  moderate  northeasterly  winds  in  the  western  half  of  the  domain  are 
advecting  somewhat  drier  air  into  the  domain,  as  can  be  seen  from  Figure  8- 
21c.  Figures  8-22  and  8-23  show  the  model  evolution  (large  scale  and  nest)  of 
low  level  temperature,  dew  point,  and  winds  for  Greenwo^  MS  (GWO)  and 
Montgomery,  AL  (MGM).  The  diurnal  temperature  evolutions  are  similar  at 
the  two  stations,  but  the  dew  point  at  MGM  drops  throughout  the  day  from  the 
dry  advecdon,  while  the  dew  point  at  GWO  rises  for  most  of  the  day. 

(2)  Both  simulations  have  a  spurious  drop  in  dew  point  in  the  afternoon,  a 
characteristic  problem  of  the  PBL  scheme.  What  seems  to  happen  is  that 
strong  mixing  through  a  significant  depth  of  the  lower  atmosphere  occurs 
suddenly  when  the  PBL  has  destabilized  enough  to  enter  the  free  convection 
regime.  When  air  which  is  relatively  dry  overlays  the  surface,  it  is  mixed  down 
to  the  surface  too  rapidly.  This  is  a  consequence  of  the  plume  model  of  the 
unstable  PBL,  which  allows  sinking  and  mixing  of  air  at  the  top  of  the  PBL  to 
the  surface  immediately.  One  possible  solution  is  to  make  a  mo^ficadon  to  the 
scheme,  which  would  allow  si^ace  parcels  to  ascend  and  mix  to  the  top  of  the 
PBL  immediately  (rising  plumes),  while  constraining  downward  mixing  to  a 
lesser  rate,  on  the  assumption  that  negatively  buoyant  parcels  do  not  traverse 
through  the  entire  depth  of  the  PBL  in  the  same  way  that  rising  parcels  do. 

(3)  The  stronger  easterly  winds  in  the  Gulf  of  Mexico  are  coming  onshore  in 
Louisiana,  leading  to  some  weak  convergence  and  light  rainshowers  through 
the  afternoon.  Ol^rvations  suggest  that  ^is  did  occur,  with  towering  cumulus 
and  occasional  showers  observ^  at  stations  in  and  near  New  Orleans. 
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Figure  8-21  MASS  model  fields  at  1800  UTC 14  September:  (a)  Temperature  ( T)  and  (b)  wit^ 
in  the  lowest  model  layer,  (c)  the  average  relative  humidity  in  the  layer  from  the  surface  to  500 
mb  (%),  and  (d)  model  precipitation  (hmdredths  of  an  inch)  accumulated  in  the  previous  three 
hours . 
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time  (UTC) 


Figure  8-22  The  diurnal  evolution  of:  (a)  temperature  (T).  (b)  dew  point  (T),  and  (c)  wind 
speed  (m  S'^)  during  14  September  1992.  at  Greenwood,  MS.  The  first  curves  (open  boxes) 
represem  the  observed  fields,  and  the  other  two  lines  are  taken  from  model  grid  points  nearest  to 
the  station  location  for  the  large  scale  (50  km)  and  nested  (15  km)  simulations  for  that  day. 


SBIR  Phase  n  Hnal  Report 


DAAD07-90-C-0134 


Page  96 


tiine(UTC) 


ws 

WS  large 
WSrresi 


Figure  S-23  The  diurnal  evolution  of:  (a)  temperature  ( T),  (b)  dew  point  (  T),  and  (c)  wind 
speed  (m  during  14  September  1992.  at  Montgomery,  /4L.  The  first  curves  (open  boxes) 
represent  the  observed fields,  and  the  other  two  lines  are  taken  from  model  grid  points  nearest  to 
the  station  location  for  the  large  scale  (50  km)  and  nested  (15  km)  simulations  for  that  day. 
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8.3.3  August  22  Southwest  Case:  Strong  Terrain  Influences 

On  the  meso-a  scale,  tenain  tends  to  be  fairly  poorly  resolved,  and  only  large  terrain 
features  which  are  well-resolved  by  the  model  grid  have  a  significant  influence  on  simulations. 
At  the  lower  limit  of  hydrostatic  mesoscale  simulations  (10-15  km),  smaller  terrain  features  and 
stronger  slopes  are  resolved,  which  have  a  great  effect  on  many  nested  simulations.  For  the  first 
set  of  dates  (August  22-24),  the  nested  New  Mexico  simuladons  used  a  10  km  resolution,  rather 
than  the  15  km  resolution  of  the  later  simulations.  The  strong  effects  of  relatively  fine-scale 
mountain  ranges  are  readily  seen.  Figure  8-24  shows  the  terrain  heights  for  this  grid.  Figures  8- 
25  through  8-27  show  a  six  hour  period  in  the  middle  of  the  nes^  10  km  simulation  for  22 
August  1992.  The  moisture  from  Hurricane  Lester  was  beginning  ro  move  in  from  the  uvst,  and 
precipitation  develops  during  die  day.  The  precipitation  forms  over  the  various  terrain  features 
of  the  area  preferentially,  fiiat  over  the  high  terrain  in  the  western  pan  of  the  grid  and  then  over 
the  Sacramento  Mountains  toward  the  east.  Temperature  perturbations  (low  temperatures  over 
the  mountains)  may  be  exaggerated  by  the  difficulty  of  formulating  a  lana^  diffusion  scheme  for 
sigma  layers  over  strongly  sloping  terrain.  Venical  (physical)  mixing  processes  typically 
b^me  entangled  with  die  horizontal  diffusion,  which  is  for  purely  numerical  purposes.  The 
horizontal  diffusion  is  designed  to  prevent  the  accumulation  of  energy  into  shon  length  scales, 
and  subsequent  aliasing  to  longer  lengths.  The  proper  treatment  of  diffusion  of  moisture  along 
sloping  surfaces  is  also  poorly  understood. 
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Figure  8-2S  Nested  model  fields  for  1500  UTC  22  August  1992:  (a)  low  level  dew  point  ( T),  (b) 
3  hr  accumulated  precipitation  (hundredths  of  an  inch),  (c)  low  level  winds  (knots),  and  (d)  low 
level  temperature  (T). 
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Figure  8-26  Nested  model  fields  for  1800  UTC  22  August  1992:  (a)  low  level  dewpoint  (  T),  (b) 
3  hr  accumulated  precipitation  (hundredths  of  an  inch),  (c)  low  level  winds  (knots),  and  (d)  low 
level  temperature  ( T). 
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Figure  8-27  Nested  model  fields  for  2100  UTC  22  August  1992:  (a)  low  level  dew  point  (  T),  (b) 
3  hr  accumulated  precipitation  (hundredths  of  an  inch),  (c)  low  level  winds  (knots),  and  (d)  low 
level  temperature  (  T). 
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8.3.4  December  10  Great  Plains  Case:  Receding  Cloud  Boundary 

The  explicit  and  implicit  treatment  of  clouds  in  mesoscale  models  remains  one  of  the 
most  difficult  scientific  and  practical  problems.  In  some  simulations,  poor  simulation  of 
cloudiness  leads  to  poor  results,  as  temperatures  are  affected  and  inland  heating  gradients  due  to 
the  delineation  of  cloudy  and  cloud-fiee  air  either  develop  when  they  shouldn't  or  fail  to  develop 
when  they  should. 

Figure  8-28  shows  the  observed  and  model  evolution  of  surface  variables  at  Columbia, 
MO  (COU)  on  10  December  1992.  In  this  case,  a  cloudy  area  associated  with  the  large  scale 
system  which  produced  the  Nonheastem  snowstorm  (Section  8.3.1)  was  receding  to  the  east. 
*nie  fairly  coti^lex  and  subtle  variations  during  the  day  were  well-simulated  by  toth  the  large 
scale  and  nested  simulations,  although  the  temperatures  were  systematically  several  degrees  too 
warm.  Figure  8-29  shows  the  modeled  clouds  (as  indicated  by  the  shortwave  transmissivity 
field,  where  low  transmissivity  indicates  thick  clouds)  moving  toward  the  northeast  in  the  nested 
IS  km  Illinois  grid.  Table  8-3  lists  the  Columbia  observations  through  the  day. 

Table  8~3  Hourly  observations  at  Columbia,  MO  for  10  December  1992. 


time 

(UTC) 

temperature  / 
dew  point 
(*F) 

wind 

speed/wind 

direction 

(deg/kts) 

clouds/sky 

conditions 

0000 

32/32 

IHIHEQSlHHi 

overcast/snow 

0100 

32/32 

140/3 

ovetcast/snow 

0200 

32/32 

BHHUESSHilH 

overcast/foe 

0300 

32/32 

180/3 

overcast/foK 

32/32 

200/4 

overcast/fOK 

i  0500 

32/32 

33/33 

200/3 

1  0700 

33/33 

220/5 

obscured/foe 

34/34 

obscured/foe 

35/35 

obscured/foe 

1  1100 

37/36 

260/7 

overcast 

« 

230/3 

clear 

1300 

36/33 

300/9 

overcast 

1400 

35/33 

310/6 

overcast/foe 

35/33 

overcast/drizzle 

1  1600  { 

35/33 

IHHlSSiHHI 

overcast/dhzzle 

36/33 

overcast/fojj 

1  1800  { 

38/34 

overcast 

HllBBiM 

40/36 

310/7 

broken 

2000 

43/35 

320/7 

scauered 

2100 

45/34 

290/7 

clear 

44/32 

280/7 

clear 

1  230ai!-"-:---:H 

42/31 

mmKMmm 

clear 

The  inferred  sequence  of  events  at  Columbia  is  as  follows; 

(1)  Conditions  slowly  improve  in  the  morning  under  southerly  or  southwesterly 
flow.  The  temperature  and  dew  points  rise. 
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(2)  A  weak  frontal  passage  occurs  at  about  1300  UTC,  with  a  wind  shift  to  the 
northwest,  and  some  dnzzle  along  the  frontal  boundary.  The  wind  speed  jumps 
up  significantly.  Ute  temperature  and  dew  points  drop  in  the  new  air  mass 
between  1200  and  1500 ITTC. 

(3)  As  skies  be^n  to  clear  out  behind  the  front  from  1500  to  2100  UTC, 
temperatures  increase  again,  and  the  dew  points  increase  also,  possibly  from 
some  evaporation  of  the  abundant  surface  moisture,  driven  by  sunshine  and 
warming. 

(4)  Drier  air  advects  in  behind  the  front,  dropping  the  dew  points  after  2000  UTC. 
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Figure  8’28  The  diurnal  evolution  of:  (a)  temperature  ( T),  (b)  dew  point  (  T),  and  (c)  wind 
speed  (m  s'^)  during  10  December  1992.  at  Columbia.  MO.  The  first  curves  (open  boxes) 
represent  the  observed  fields,  and  the  other  two  lines  are  taken  from  model  grid  points  nearest  to 
the  station  location  for  the  large  scale  (50  km)  and  nested  (15  km)  simuUuions  for  that  day. 
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Figure  8-29  The  MASS  model  shortwave  transmssivity  field  at:  (a)  1500  UTC  (b)  1800  UTC, 
and  (c)  2100  UTC  10  December  1992.  Low  values  of  transmissivity  indicate  the  presence  of 
optically  thick  cloud  cover,  either  explicitly  simulated  or  irf erred  from  the  relative  humidity 
field. 
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9.  Development  of  Graphical  Interface 

The  graphical  interface  for  the  MASS  simulation  system  was  originally  planned  to  be 
designed  for  a  NeXT  Color  Cube  system  using  an  Interface  Builder  under  the  NeXTStep 
environment.  Because  of  the  use  of  Motif  for  the  IMETS  system  at  the  ASL,  there  was  a  shift  in 
the  development  of  the  graphical  interface  from  the  NeXTStep  environment  to  the  X 
Windows/Motif  environment.  This  approach  has  produced  an  X  Windows/Motif-based  interface 
that  will  be  consistent  with  the  Army's  software  standards,  and  hopefully,  one  that  can  be 
implemented  as  part  of  the  IMETS  system  in  the  funire  if  desired.  This  development  was  done 
on  the  Stardent  7S0  which  was  delivered  to  ASL. 


9.1  Design  of  the  GUI 

The  Graphical  User  Interface  (GUI)  was  developed  using  the  Motif  Resource  Manager 
(MRM),  which  allows  widgets  (i.e.  bunon  and  text  boxes)  to  be  created  based  on  information 
contain^  in  separate  input  files  called  User  Interface  Definition  (UID)  files.  The  information 
contained  in  these  files  includes  the  location  and  sizes  of  different  features,  text  font  choices, 
labels  and  titles  and  default  initial  choices.  Ihese  files  are  separate  from  the  main  source  code 
written  in  C,  which  compiles  very  slowly  due  in  part  to  the  large  X  and  Motif  libraries  needed. 
The  separation  of  the  UID  files  from  the  C  source  code  allows  the  developer  to  dramatically 
improve  his  efficiency  by  modifying  the  UID  files  rather  than  the  C  files  when  designing  the 
structure  of  the  windows.  Features  can  quickly  be  re-sized  and  moved  around  or  renarr^  in  the 
UID  files  and  tested. 

The  GUI  was  designed  so  that  the  user  can  modify  and  piod,  '«  the  necessary  input  files 
for  the  preprocessor  and  model  as  well  as  submit  preprocessor,  nudging  or  model  runs 
completely  from  the  screen.  The  design  of  the  menus  is  shown  in  Figure  9-1.  The  top  menu 
(Figure  9-2)  has  options  for  configuring  either  the  preprocessor  or  the  mass  option  files  or  setting 
up  and  submitting  a  simulation.  If  the  user  chooses  to  configure  input  option  files,  subsequent 
menus  let  the  user  choose  which  specific  option  file  he  would  like  to  generate.  These  option  files 
will  be  created  and  stored  in  the  current  directory  of  the  user  at  the  time  the  GUI  was  staned. 
There  are  six  configuration  menus  for  the  preprocessor  as  listed  in  Figure  9- 3a,  for  the  six 
different  inodules:  prepgrd,  prepdat,  prepro,  prepbog,  prepbc  and  prepnudg.  There  are  three 
configuration  menus  for  the  MASS  model  (Figure  9-3b),  one  for  the  MASS  model  and  two  for 
the  SRPH  and  nudging  sections  of  the  MASS  option  file.  The  configuration  of  any  of  these  three 
menus  results  in  the  generation  of  the  same  mass. opt  option  file.  If  the  user  chooses  to  configure 
and  submit  a  simulation,  the  submit  menu  (Figure  9-4)  allows  the  user  to  do  any  or  all  of  the 
following: 

(1)  Choose  the  Preprocessor  Init  Tag  and  run  any  or  all  of  the  preprocessor 
modules. 

(2)  Run  the  nudging  module  for  up  to  six  different  nudging  times. 

(3)  Choose  the  Model  Run  Tag  and  run  the  MASS  model. 

The  configuration  menus  for  the  six  different  preprocessor  modules  are  shown  in  Figures 
9-5  through  9-10.  These  menus  allow  the  user  to  completely  define  the  input  parameters  needed 
by  the  preprocessor  to  run  a  simulation.  The  three  configuration  menus  needed  to  run  a  MASS 
simulation  are  shown  in  Bgures  9-11  through  9-13. 
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Figure  9~l  Schematic  diagram  of  the  basic  graphical  interface  structure. 
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Figure  9-2  The  mem  in  the  imass  interface. 
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Figure  9-3  The  top  menu  lists  for  the  :  (a)  preprocessor  and  (b)  MASS  model. 
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PREPGRD  Module  Conllguntion  Menu 
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Figure  9~5.  The  prepgrd  cor^guration  menu. 
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Figure  9-1.  The  prepro  configurcaion  menu. 
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Figure  9S.  The  prepbog  cor^guration  menu. 
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Figure  9-9.  The  prepbc  corrfiguration  menu. 
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Figure  9-10.  The  prepnudg  configuration  menu. 
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SRPH  Module  Conf^uration  Menu 
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Figure  9~I2.  The  SRPH  coi^guration  menu. 
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Figure  9-13.  The  Nudging  corfiguration  menu. 
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10.  Summary  and  Conclusions 

The  fundamental  achievement  of  this  SBIR  project  was  the  creation  of  a  self-contained 
mesoscale  atmospheric  simulation  system  which  can  generate  real-time  or  historical  simulations 
on  a  high  perfOTinance  moderate-cost  workstadon  computer.  The  system  can;  (1)  ingest  a  variety 
of  atmospheric  data  types;  (2)  combine  the  diverse  mixture  of  ingested  data  into  a  dynamicdly 
consistent  initializadon  dataset  for  the  numerical  model;  (3)  generate  a  mesoscale  numerical 
simuladon;  and  (4)  display  the  ouq)ut  from  the  simuladons  in  a  variety  of  formats.  A  review  of 
the  signiHcant  accomplishments  during  the  development  and  evaluation  of  the  system  is 
presented  in  Secdon  10.1.  A  discussion  of  the  areas  in  which  addidonal  research  is  needed  to 
improved  the  mesos^e  simuladon  system  is  presented  in  Secdon  10.2. 


10,1  Project  Summary 

The  project  ccmsisted  of  three  significant  components.  In  the  first  portion  of  the  project, 
the  workstadon-t^ed  mesoscale  atmospheric  simuladon  system  was  created  by  porting  a  version 
of  the  Mesoscale  Atmospheric  Simulation  System  (MASS)  from  a  supercomputer  to  a  Stardent 
750  vector  processing  workstation,  and  then  upgrading  the  system  by  implementing:  (1)  a  new 
lateral  boundary  condition  scheme;  (2)  a  positive  definite  ^vection  formulation;  (3)  a  more 
sophisticated  surface  energy  and  moisture  budget  formulation;  (4)  a  four-dimensional  data 
assimilation  system  based  on  a  Newtonian  relaxation  scheme;  and  (S)  a  method  to  enhance  the 
initialization  of  relative  humidity  through  the  use  of  satellite  image  data,  MDR  reports,  pilot 
reports  and  surface  cloud  observations.  In  the  second  portion  of  the  project,  a  graphical  user 
interface  (GUI)  was  developed  to  permit  the  user  to  easily  reconfigure  the  system  and  execute 
simulations,  and  the  model  software  was  modified  to  increase  its  computational  performance  on 
the  workstation  computer.  The  final  segment  of  the  project  consisted  of  the  execution  of  a 
observation  system  simulation  experiment  (OSSE)  to  test  the  performance  of  the  data 
assimilation  system,  and  the  execution  of  36  real  data  simulations  to  evaluate  the  performance  of 
the  simulation  system  in  a  variety  of  environments. 

The  first  improvement  to  the  modeling  system  in  this  project  was  the  implementation  of  a 
radiative  lateral  boundary  condition  formulation.  The  formulation  was  based  on  the  scheme 
described  by  Orlanski  (1976).  In  this  scheme,  atmospheric  features  which  approach  the  lateral 
boundaries  are  permitted  to  propagate  through  the  boundary  by  calculating  a  composite  phase 
velocity  at  each  boundary  point  and  allowing  the  features  to  propagate  through  the  lateral 
bounds^  with  this  phase  velocity.  This  foimulation  improved  the  qu^ity  of  the  simulations  near 
the  latei^  boundaries  and  also  rkluced  the  portion  of  the  nxxlel  domain  (i.e.  the  number  of  grid 
points)  devoted  to  the  implementation  of  the  boundary  conditions.  Under  the  previous 
Kreitzberg-Peikey  sponge  boundary  condition,  four  boundary  rows  and  columns  were  required  to 
blend  the  external  and  internal  tendencies  and  filter  the  outwardly  propagating  features  from  the 
model  domain.  This  effectively  reduced  the  useful  portion  of  the  model  matrix  (i.e.  the  portion 
in  which  the  evolution  of  the  prognostic  variables  is  determined  solely  by  the  model  physics)  by 
8  grid  points  along  each  coordinate  axis.  In  contrast,  the  radiative  scheme  utilizes  only  the 
outermost  row  and  column  for  the  specification  of  the  boundary  conditions  and  there  is  no 
necessity  to  have  a  highly  diffusive  sponge  region  for  several  rows  or  columns  adjacent  to  each 
lateral  boundary. 

The  quality  of  the  modeling  system  was  also  improved  in  this  project  by  the 
implementation  of  a  positive  definite  honzontal  advection  scheme  with  low  implicit  diffusion  for 
the  advection  of  liquid  and  frozen  water  substances  and  any  passive  tracer  substances  that  the 
user  decides  to  include  in  a  model  simulation.  The  advection  scheme  is  a  version  of  the 
MPDAT^A  scheme  originally  formulated  by  Smolaikiewicz  (1983a,  1983b).  This  scheme  has  the 
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ability  to  maintain  shaxp  boundaries  in  simulated  fields  while  not  generating  furious  numerical 
oscillations  which  cause  unrealistic  negative  values  to  appear  in  proximity  to  regions  of  strong 
gradients.  The  scheme  was  extensively  tested  in  idealized  one-diiMiisional  and  two-dimensional 
simulations  before  it  was  implemented  into  the  three-dimensional  version  of  the  MASS  model. 
It  will  be  possiUe  u>  implement  this  advection  scheme  for  all  of  the  iiKxiel  prognostic  variables 
in  the  future  with  a  modest  additional  effort 

Various  pans  of  the  sioface  and  planetary  boundary  layer  parameterization  schemes  were 
signiHcantly  improved  as  a  part  of  this  project.  A  better  method  of  calculating 
evaptmanspiration  firom  land  surfaces  was  formulate  in  which  the  transpiration  frmn  plants  is 
sensitive  to  the  fractional  vegetation  cover  (inferred  from  renootely-sensed  vegetation  index  data) 
and  land  use  type,  as  well  as  soil  moisture.  Evaporation  is  also  parameterized  from  tnue  soil,  and 
from  a  cover  moisture  reservoir  consisting  of  intercepted  rainfall  and  dewfall.  As  a  pan  of  the 
evapotran^iration  effort,  the  model  hydrology  framework  was  extended  to  two  soil  layers  and 
the  cover  reservoir,  and  paiameterizations  for  rainfall  interception  and  snow  cover  were  atkied. 

Tests  were  made  with  the  initialization  of  surface  and  subsoil  temperature.  It  was  found 
that  the  model  is  quite  sensitive  to  subsoil  temperanire,  especially  in  the  simulation  of  nighttime 
low  temperatures.  Using  an  average  temperature  over  the  previous  few  days  seems  to  be  a 
reasonable  way  to  initialize  the  subsoil  temperature. 

A  formulation  for  the  inclusion  of  slope  effects  on  the  shortwave  radiation  received  at  the 
surface  was  also  added  to  the  model.  The  local  terrain  slope  magnitude  and  slope  azimuth  are 
used  to  calculate  a  corrected  solar  zenith  angle,  which  replaces  the  actual  zenith  angle  in  the 
radiation  equations.  For  high  resolution  nest^  simulations,  the  grid-resolved  terrain  slopes  can 
be  large  enough  to  significantly  alter  the  local  radiation  budget. 

The  initialization  of  the  mesoscale  simulation  model  was  improved  in  this  project  by 
impletnenting  four  types  of  four  dimensional  dau  assimilation  (FDDA)  schemes  and  an 
innovative  method  of  static  moisture  initialization.  These  FDDA  schemes  include  (1)  Newtonian 
relaxation  or  nudging  of  gridded  rawinsonde,  surface  or  profiler  observations  (Hoke  and  Anthcs, 
1976)  in  which  the  tnodel  state  is  relaxed  towards  an  analysis  of  observations;  and  (2)  the  use  of 
Manually  Digitized  Radar  (MDR)  data  to  specify  moisture  convergence  (precipitation  rates)  in 
areas  which  are  (are  not)  subject  to  convection  according  to  the  Kuo-MESO  cumulus 
parameterization  scheme,  llie  static  moisture  initialization  is  based  upon  the  enhancement  of  the 
three  dimensional  moisture  analysis  through  the  use  of  surface  observations  of  clouds,  pilot 
reports,  MDR  data  and  infrared  satellite  images.  The  surface,  rawinsonde  and  profiler  nudging 
schemes  are  similar  to  the  schemes  used  by  Staufrer,  et  al.,  (1991).  The  most  significant  changes 
from  Stauffer,  et  al.,  (1991)  are  an  improvement  in  the  method  of  calculating  the  analysis 
confidence  factors  and  a  scheme  to  improve  the  surface  analysis  near  coastlines.  Temperatures 
are  derived  from  the  profiler  wind  field  through  the  inversion  of  the  divergence  equation  (Cram, 
et  al.,  1991)  and  assimilated  along  with  the  wind  field. 

Simulation  experiments  indicated  that  the  assimilation  of  MDR  data  helped  to  accurately 
define  both  the  grid  scale  and  subgrid  scale  precipitation  field  at  the  time  of  model  initialization. 
The  benefits  of  MDR  assimilation  typically  lasted  for  about  6  to  9  hours.  They  were  most 
dramatic  during  the  initial  three  hours  of  a  simulation  due  to  the  elimination  of  the  precipitation 
spin-up  problem  that  is  common  to  all  numerical  models. 

The  synthetic  RH  scheme  is  most  valuable  when  simulating  cases  with  weak  or 
nonexistent  synoptic  forcing  and  strong  boundary  layer  forcing.  In  these  situations,  convection  is 

often  triggered  by  differential  heating  due  to  meso-|S  scale  variations  in  cloud  cover.  The 
synthetic  RH  scheme  is  the  component  of  the  model  data  assimilation  system  which  is  most 
effective  at  adding  such  frne  scale  variations  in  cloudiness. 

The  results  of  the  observing  system  simulation  experiment  (OSSE)  indicated  that  the 
main  benefit  of  nudging  is  the  reduction  of  the  spin-up  problem  during  the  initial  3  to  6  hours  of 
a  simulation.  After  this  time,  nudged  simulations  were  inferior  to  a  statically  initialized 
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simulation  because  the  nudged  simulation  had  to  be  initialized  12  hours  earlier.  The  poor 
peiformance  of  the  nudged  simulations  after  6  hours  was  most  likely  due  to  the  fact  that  the 
boundary  conditions  were  based  on  an  NGM  12  hour  forecast  rather  than  "observations"  from  a 
surrogate  atmosphere  simulation  (SAS).  During  the  assimilation  period  and  the  first  6  hours  of 
the  forecast,  this  infericMr  boundary  condition  information  propagated  across  a  large  enough 
portion  of  the  model  ^main  to  (tegrade  the  simulation.  The  results  of  the  case  of  5-6  Septemter 
1992  indicate  that  the  assimilation  of  rawinsonde  data  produces  the  ^eatest  improvement,  while 
surface  data  is  less  valuable  and  profiler  data  suffers  from  a  lack  of  information  abcnit  the  mass 
field.  Since  the  OSSE  consisted  of  only  1  case,  these  results  are  not  necessarily  true  for  all 
conceivable  situations.  For  cases  with  weak  forcing,  the  assimilation  of  rawinsonde  data  may  be 
less  important  while  surface  data  assimilation  may  have  a  greater  impact.  Stauffer  et  al.  (1^1) 
found  tins  to  be  true  in  their  weak  synoptic  forcing  case. 

Given  an  abundance  of  time  and  other  resources,  OSSEs  have  the  potential  to  expand  the 
ut^erstanding  of  data  assimilation  beyond  the  limited  results  presented  here.  Several  potentially 
viduable  experiments  must  wait  for  another  opportunity.  They  include: 

1)  the  use  of  observations  from  the  SAS,  rather  than  output  from  a  numerical 
model,  to  provide  boundary  condition  information.  This  may  help  to  reduce  the 
degradation  of  the  simulation  which  occurs  when  forecast  boundary  condition 
dau  propagates  into  the  domain  during  the  data  assimilation  period. 

2)  creating  a  technique  to  simulate  MDR,  IR  satellite  and  other  non-standard  data 
from  the  SAS  so  that  the  MDR  assimilation  scheme  and  the  synthetic  RH 
scheme  could  be  tested. 

3)  testing  the  effects  of  the  temporal  and  spatial  distribution  of  the  simulated 
observations  on  the  numerical  forecast. 

Quite  a  large  number  of  MASS  simulations  were  performed  for  this  project,  including 
both  the  OSSE  simulations  and  those  simulations  performed  as  part  of  the  36-case  sample.  Table 
10-1  summarizes  the  types  of  phenomena  that  the  model  tended  to  simulate  well  and  those  types 
of  features  which  were  consistently  difficult  to  simulate.  Some  of  the  strengths  and  weaknesses 
are  characteristic  of  mesoscale  modeling  in  general,  and  some  are  peculiar  to  MASS.  Model 
problems  can  be  separated  into  two  main  categories:  (1)  those  which  arise  from  a  lack  of  initial 
or  boundary  data  of  sufficient  resolution,  frequency,  or  quality;  and  (2)  those  which  are  caused 
by  shortcomings  in  the  model  physics,  arising  either  from  inadequate  understanding  of  the 
physical  process,  or  from  constraints  imposed  by  the  model  numerical  formulation  or  by  model 
chaiacteristics  such  as  the  grid  spacing.  Convective  and  cloud  phenomena  are  the  areas  in  which 
mesoscale  models  can  make  the  greatest  improvement  over  current  operational  models.  Despite 
some  successes  in  these  areas,  profound  difficulties  remain  in  the  parameteiization  of  convection 
and  other  cloud  processes. 
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Table  10-1  Types  cf  meteorological  phenomena  which  were  generally  well-simulated  or  difficult 
to  simulate  with  the  Mesoscale  Atmospheric  Simulation  System. 


M  1  1  III  rill  f'lTT"!  1  T1 

•Synoptic  or  meso-a  scale  features 
which  are  well-resolved  by  the 
lawinsoode  networic. 

•Convective  events  which  are  only 
weakly  forced  by  larger  scale 
circulations,  i.e.  quasi-barotropic 
systems. 

•Convective  events  which  are  strongly 
forced  by  larger  scale  circulatioas. 

•Features  which  result  from  convective 
feedback  to  the  grid  scale.  e.g. 
from  latent  heating,  downdrafts, 
cool  outflows,  etc. 

•Mesoscale  circulations  which  are  tied 
to  surface  features  that  are  well- 
rqxesented.  e.g.  land/sea  breezes, 
terrain-induced  chculations. 

•  Mesoscale  circulations  which  are 
driven  by  more  subtle  surface 
gradients,  eqKdally  soil  moisture. 

•  Surface  and  boundary  layer  flows 
occurring  under  quiescent  (clear, 
non-convective)  conditions. 

•Circulations  driven  by  cloud  boundary 
effects. 

•Features  in  the  vicinity  of  large  terrain 
gradients. 

An  alwayS'important  issue  in  numerical  weather  prediction  is  the  computational  time 
required  to  make  simulations.  With  the  many  different  model  options,  a  wide  range  of  speeds 
can  be  obtained  on  a  given  grid,  depending  on  the  model  physics  selected,  and  the  frequency  at 
which  parameterization  schemes  are  invoke.  Hie  Stardent  750  Workstation  has  proven  to  be  an 
effective  platform  for  the  MASS  preprocessor  and  model,  and  the  ability  to  add  another 
processor  would  increase  the  speed  considerably.  At  the  end  of  the  36-case  sample  simulations, 

the  time  requited  to  make  a  24  hr  large  scale  run  on  a  55  x  50  domain  with  a  grid  spacing  of  SO 
km  was  about  3.7  hr.  The  nested  simulations  take  longer,  because  of  the  much  shorter  timestep: 

20  s  versus  75  s  for  the  large  scale  runs.  A  12  hr  nested  simulation  on  a  55  x  50, 15  km  grid  took 
about  4.9  hr  on  the  Stardent  750.  These  times  are  a  little  slow  for  operational  use,  but  the 
Stardent  line  of  workstations  is  more  than  two  years  old,  and  much  faster  machines  are  now 
available.  A  great  advantage  of  the  Stardent  is  the  ability  to  perform  vector  processing,  a 
capability  which  is  currently  available  only  on  much  more  expensive  computers,  and  is  not 
available  on  fast,  new  workstations  such  as  the  Hewlett-Packard  9000  series  (the  "Snake")  and 
the  new  "Alpha"  woikstations  from  Digital  Equipment  Corp.  These  machines  are  much  faster  on 
scalar  calculations  than  the  Stardent's  MIPS  R3000  processor,  but  it  remains  to  be  seen  how 
much  is  lost  by  the  absence  of  vector  processing.  Some  parts  of  the  MASS  model  (advection, 
diffusion)  vectorize  very  well,  significantly  speeding  up  the  execution  of  the  code.  On  the  other 
hand,  some  parts  of  the  model  (surface  and  microphysical  parameterizations)  do  not  vectorize 
well,  and  would  benefit  greatly  from  faster  scalar  processing  power.  Figure  10-1  shows  the 
percentage  of  execution  time  required  for  the  major  parts  of  the  model  when  the  code  is  compiled 
with  the  vectorizadon  option.  The  most  time-consuming  single  component  of  the  model  is  the 
SRPH  scheme  (Surface  energy  budget-Radiauon-Planetary  boundary  layer-Hydrology). 
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Figure  10^1  Percentage  of  model  execution  time  required  for  major  parts  of  the  MASS  model. 
The  run  used  for  this  breakdown  was  a  45  bn  simulation  on  the  Stardent  750,  with  nudging  off, 
the  prognostic  microphysical  scheme  on,  and  the  Kuo-MESO  convective  parameterization 
scheme  on.  All  the  modules  were  optinUzed  with  the  vectorization  option  of  the  Stardent 
FORTRAN  compiler. 

102  Areas  for  Future  Development 

Despite  the  substantial  progress  that  has  been  made  in  the  develt^ment  of  mesoscale 
models  during  the  past  decade,  and  the  widespread  use  of  these  systems  in  recent  years,  there  are 
still  a  numbCT  of  problem  areas  diat  need  to  be  addressed  before  mesoscale  simulations  can 
achieve  their  ultimate  potential.  The  ma^  issues  that  n^d  to  be  addressed  include:  (1)  the 
parameterization  of  moist  convection  and  its  interaction  with  grid  scale  moisture  physics;  (2)  the 
parameterization  of  the  boundary  layer,  (3)  the  modeling  of  surface  processes;  (4)  the 
representation  of  clouds  and  their  interaction  with  radiative  processes;  (S)  the  availability  of 
simicient  data  to  define  mesoscale  features  in  the  initial  state;  (6)  the  assimilation  of  data  tom 
new  observing  systems  into  the  mesoscale  model  in  a  beneficial  manner,  and  (7)  the  endless  need 
for  higher  computational  performance.  A  brief  review  of  each  one  of  these  problem  areas  will 
provide  an  indication  of  the  direction  of  the  future  develq[>ment  of  MASS  and  other  mesoscale 
models 
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One  of  the  most  critical  problems  is  the  parameterization  of  the  grid  scale  effects  of  the 
sub-grid  scale  processes  associated  with  shallow  and  deep  moist  convectitm  in  mesoscale  models 
with  g^  increments  tetween  approximately  S  km  and  20  km.  A  review  of  the  inany  unresolved 
issues  of  convective  parameterization  at  this  scale  has  been  compiled  by  Molinari  and  Dudek 
(1992).  At  this  resolution,  the  mesoscale  model  partially  resolves  the  convecdve-scale 
circulation  but  caiuiot  resolve  the  intense  updrafts  and  downdraft  which  accomplish  much  of  the 
vertical  transport  of  energy,  moisture  and  momentum.  Tbus,  the  parameterization  scheme  must 
simulate  the  gr^  scale  ejects  of  the  large  vertical  transports  which  occur  in  the  updrafts  and 
dowiKlrafts,  but  not  duplicate  the  effects  of  the  processes  which  are  resolved  on  the  mesoscale 
grid.  If  a  parameterization  is  not  used  at  this  scale,  the  absence  of  the  vertical  transport  by 
convective-scale  uptfarafts  and  downdrafts  fr^uently  results  in  the  spurious  feedback  process 
which  can  unrealistically  amplify  convective  systems.  The  parameterized  updrafts  and 
downdrafts  work  against  the  spurious  feedback  by  transporong  heat  and  moisture  upward.  The 
proper  approach  to  the  parameterization  of  convection  at  this  scale  is  still  subject  to  considerable 
debate.  Most  mesoscale  nxxlelers  have  simply  utilized  schemes  such  as  those  formulated  by 
Kuo  (1965)  and  Anthes  (1977),  that  were  intended  for  use  in  noodels  with  grid  spacing  larger 
than  20  lom.  Two  schemes  were  expressly  designed  for  models  with  grid  increments  in  the  20  to 
25  bn  range:  the  schemes  of  Fritsch  and  Chappell  (1980)  and  Frank  and  Cohen  ( 1987).  A  mixed 
set  of  results  has  been  documented.  An  appuling  approach  is  to  avoid  the  entire  moist 
convection  parameterization  issue  by  explicitly  simulating  the  convection  with  a  non-hydrostatic 
cloud  scale  model  executed  over  a  mesoscale  dmnain.  A  few  experimental  simulations  of  this 
type  have  been  performed  with  a  non-hydrostatic  version  of  the  CSU  RAMS  model,  but 
computational  requirements  will  restrict  such  simulations  to  very  small  domains  for  the 
imm^ate  future.  Thus,  the  improvement  in  convective  parameterization  schemes  will  play  an 
important  role  in  determining  the  skill  of  a  mesoscale  model  to  simulate  the  evolution  of 
mesoscale  convective  systems.  However,  cloud-scale  simulations  of  mesoscale  systems  will 
undoubtedly  provide  an  important  dataset  with  which  to  understand  how  to  parameterize 
convective  systems  in  mesoscale  models. 

The  parameterization  of  the  boundary  layer  is  another  area  in  which  mesoscale  models 
can  be  improved.  A  key  component  of  the  boundary  layer  parameterization  problem  is  the 
interaction  of  the  model's  boundary  layer  parameterization  with  shallow  and  deep  convective 
clouds  that  are  rooted  within  the  boundary  layer.  However,  most  boundi^  layer 
parameterization  schemes  do  not  directly  interact  with  clouds  that  extend  into  the  middle  and 
upper  troposphere. 

The  representation  of  surface  processes  is  also  an  important  contributor  to  the  skill  of  a 
mesoscale  simulation.  A  key  issue  in  this  area  is  the  calculation  of  the  area-averaged  surface 
moisture  flux  over  land  surfaces.  The  flux  of  moisture  from  the  surface  of  the  earth  into  the 
atmosphere  is  the  result  of  direct  evaporation  from  the  surface  soil  layer  and  the  material  on  top 
of  the  soil  (e.g.  vegetation,  pavement)  as  well  as  the  transpiration  from  vegetation  which  extracts 
water  from  the  soil  layer  below  the  surface  layer.  The  time  dependent  area-averaged 
transpiration  from  a  plant  canopy  is  a  very  difficult  quantity  to  calculate  because  of  the 
heterogeneous  transpiration  rates  of  plants  resulting  from  (1)  their  different  locations  in  the 
canopy;  (2)  the  intrinsic  physiological  properties  of  different  species;  and  (3)  large  variations  in 
soil  properties  and  moisture  contents  over  small  distances  underneath  the  canopy.  Surface 
heterogeneity  also  poses  a  significant  problem  to  the  calculation  of  area-averaged  direct 
evaporation.  A  significant  obstacle  to  the  accurate  modeling  of  this  process  is  the  lack  of  quality 
observed  evaporation  and  transpiration  data  over  a  wide  variety  of  surfaces.  Field  experiments 
such  as  FIFE  (Sellers  et  al.,  1988)  have  gathered  detailed  measurements  of  evapotranspiration 
over  only  a  few  types  of  mostly  homogeneous  surfaces. 

The  simulation  of  clouds  and  their  interaction  with  shon  wave  and  long  wave  radiation 
can  be  crucial  to  the  quality  of  a  mesoscale  simulation.  For  example,  the  quality  of  a  simulation 
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of  low  level  circulations  that  develop  due  to  the  differential  heating  of  the  surface  between 
cloudy  and  clear  regions,  is  heavily  dependent  on  the  model's  ability  to  accurately  simulate 
radiative  processes  in  the  cloudy  ai^  clear  regions.  There  are  two  significant  issues  related  to 
this  problem  area.  The  first  issue  is  the  representation  of  clouds  in  the  mesoscale  model.  One 
criteria  to  determine  the  presence  of  clouds  in  a  model  layer  is  die  existence  of  significant 
anoounts  of  simuUued  cloud  water  (ice)  or  rain  water  (snow).  However,  cloud  water  (ice)  ot  rain 
water  (snow)  can  only  exist  in  a  gpd  cell  for  a  significant  amount  of  time  if  the  air  in  the  cell  is 
saturau^.  This  means  that  the  entire  grid  cell  will  either  be  cloudy  (if  the  grid  scale  is  saturated) 
OT  clear  if  this  is  the  only  mechanism  for  cloud  specification  in  the  model.  This  is  clearly  not 
realistic  at  the  mesoscale  since  a  significant  amount  of  cloudiness  can  exist  in  a  grid  cell  even 
when  the  grid  cell  average  relative  humidity  is  somewhat  below  100%.  This  is  typically 
represented  in  mesoscale  models  by  using  a  cloud  fraction-relative  humidity  relationship  to 
spwify  a  fractional  coverage  of  clou^  in  layers  where  the  relative  humidity  is  below  100%  but 
above  a  threshold  value  (e.g.  80%).  These  clouds  are  considered  to  be  sub-grid  in  scale  and 
hence  are  not  represent^  in  the  grid  scale  cloud  water  (ice)  or  rain  water  (snow)  fields. 
However,  these  clouds  do  interact  with  the  radiative  parameterization  and  thus  can  have  a 
significant  impact  on  the  surface  energy  budget.  Unfortunately,  the  relationship  between  layer 
cloud  fraction  and  relative  humidity  is  subject  to  a  considerable  amount  of  scatter  due  to  the 
multiscale  processes  which  generate  and  destroy  clouds.  It  will  most  likely  be  necessary  to 
utilize  other  parameters  in  a^idon  to  reladve  humidity  to  obtain  a  reasonable  esdmate  of  the 
fraction  of  a  model  layer  that  is  covered  by  clouds.  A  considerable  amount  of  additional  research 
needs  to  be  done  on  the  relationship  of  grid  scale  variables  to  layer  cloud  fractions.  If  this  is  not 
a  difficult  enough  problem,  it  must  be  remembered  that  in  addition  to  the  fractional  coverage, 
the  t^tical  depth  (i.e.  the  liquid  or  frozen  water  content)  of  the  parameterized  clouds  must  be 
known  in  order  to  calculate  the  radiative  effects  of  the  clouds.  A  second  issue  relates  to  the  way 
in  which  the  radiative  calculations  are  performed.  The  execution  of  a  comprehensive  radiative 
transfer  formulation  within  a  mesoscale  model  is  not  computationally  feasible  at  present  or  in  the 
near  future.  Therefore,  simplified  parameterizations  must  be  used.  The  challenge  is  to  create  a 
computationally  efficient  radiative  scheme  that  can  simulate  the  important  aspects  of  the 
radiative  effects  of  both  parameterized  and  explicit  clouds  in  a  mesoscale  model. 

A  long  standing  problem  for  mesoscale  modelers  has  been  the  lack  of  data  with  which  to 
define  mesoscale  circulations  at  the  time  of  initialization  of  a  mesoscale  model.  In  the  past, 
rawinsonde  data  has  been  the  foundation  for  the  preparation  of  the  initialization  datasets  for 
mesoscale  models.  With  an  average  separation  of  about  400  km  between  observing  sites,  the 

rawinsonde  network  cannot  resolve  meso-P  scale  features  and  can  provide  only  a  marginal 

representation  of  meso-a  scale  systems.  Surface,  satellite  and  aircraft  data  have  provided  some 
information  about  finer  scale  mesoscale  features,  but  it  is  still  true  that  if  a  feature  is  not  well- 
resolved  by  the  rawinsonde  data  it  generally  is  not  well-represented  in  the  initial  state.  The 
success  of  mesoscale  models  under  these  conditions  has  been  the  result  of  the  fact  that  a 
significant  fraction  of  mesoscale  features  result  from  (1)  the  non-linear  interactions  between 
coarse  scale  (resolvable  by  the  rawinsonde  network)  features;  and  (2)  the  forcing  supplied  by 
small  scale  features  of  the  earth's  surface  (e.g.  land/water  boundaries,  terrain  elevation  features, 
etc.)  which  are  fixed  in  time  and  can  be  accurately  mapped.  Fonunately,  new  observing  system 
technology  is  beginning  to  provide  operational  tools  which  can  be  used  to  make  routine 
measurements  of  mesoscale  features.  These  systems  offer  the  potential  to  improve  the 
initialization  of  mesoscale  models  and  to  provide  mesoscale  data  to  improve  the  parameterization 
schemes  used  in  mesoscale  tiKxlels.  Examples  of  these  new  systems  are  Doppler  wind  profilers 
and  NEXRAD  Doppler  radars.  Unfortunately,  virtually  all  of  the  new  observing  systems  provide 
high  quality  and  high  resolution  measurements  of  only  a  subset  of  the  variables  which  are 
requir^  for  initialization.  For  example,  the  NEXRAD  system  provides  only  reflectivity  and 
radial  wind  data.  However,  a  mesoscale  model  needs  both  components  of  the  horizontal  wind 
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field,  the  temperature,  surface  messuie,  water  vapor  mixing  ratio  and  the  cloud  water  and  rain 
water  mixing  ratio  values.  The  challenge  during  the  coming  years  will  be  to  create  an 
initializatitMi  procedure  that  utilizes  a  corobinadmi  of  remote  sensing  systems  and  innovative  data 
retrieval  techniques  to  buUd  an  initialization  dataset  which  can  resolve  the  three-dimensional 

structure  of  meso-^  scale  features  such  as  mesoscale  convective  systems. 

Another  issue  that  has  frustrated  all  atmospheric  modelers  is  the  insatiable  thirst  that 
atmospheric  models  have  for  computational  power.  As  computers  become  more  powerful, 
atmosphoic  models  have  quickly  consumed  all  of  the  increased  power  by  using  higher  resolution 
grids  and  incorporating  more  detailed  physics.  Thus,  the  processing  time  for  the  "best" 
simulation  always  tends  to  be  the  maximum  acceptable  time  for  a  particular  research  or 
cqperational  iq>plication.  Fortunately,  coii:q>uter  technology  is  advancing  rapidly  at  the  present 
time.  It  is  becoming  increasingly  evident  that  the  principal  scientific  computing  engines  of  the 
mid  to  late  1990's  will  be  chancterized  by  massive  parallelism  (hund^s  to  thousands  of 
processes)  in  a  distributed  memory  configuration.  An  overview  of  the  cunent  state  of  parallel 
computing  and  prospects  for  the  future  is  presented  by  Poutain  and  Bryan  (1992).  All  the 
evidence  indicates  that  computational  performance  will  increase  dramadcally  in  all  price  ranges 
over  the  next  several  years.  The  increased  performance  will  provide  an  t^rpor^ity  for  the  user 
with  modest  resources  to  run  a  three-dimensional  mesoscale  simulation  over  a  signiricant  domain 
in  less  than  an  hour  on  a  low-cost  desktt^  computer  capable  of  generating  processing  speeds  of 
1(X)  Megaflops.  On  the  other  end  of  the  specttum,  the  user  with  substantial  resources  will  be 
able  to  execute  non-hydrostadc  cloud-scale  resolution  simulations  over  mesoscale  domains  in 
near  real-tiine  on  high  performance  massively  parallel  supercomputers  that  can  attain  processing 
speeds  between  SOO  Gigaflops  and  1  Teraflop. 
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